Kernel Support Vector Machines for Classification and Regression in C#

April 27, 2010 / cesarsouza / 53 Comments

Kernel methods in general have gained increased attention in recent years, partly due to the grown of popularity of the Support Vector Machines. Support Vector Machines are linear classifiers and regressors that, through the Kernel trick, operate in reproducing Kernel Hilbert spaces and are thus able to perform non-linear classification and regression in their input space.

Download the classification sample application (with source code)
Download the regression sample application (with source code)
Browse the framework’s examples (includes linear and logistic machines)

Foreword

If you would like to use SVMs in your .NET applications, download the Accord.NET Framework through NuGet. Afterwards, creating support vector machines for binary and multi-class problems with a variety of kernels becomes very easy. The Accord.NET Framework is a LGPL framework which can be used freely in commercial, closed-source, open-source or free applications. This article explains a bit how the SVM algorithms and the overall SVM module was designed before being added as part of the framework.

Introduction

Support Vector Machines

Support vector machines (SVMs) are a set of related supervised learning methods used for classification and regression. In simple words, given a set of training examples, each marked as belonging to one of two categories, a SVM training algorithm builds a model that predicts whether a new example falls into one category or the other. Intuitively, an SVM model is a representation of the examples as points in space, mapped so that the examples of the separate categories are divided by a clear gap that is as wide as possible. New examples are then mapped into that same space and predicted to belong to a category based on which side of the gap they fall on.

A linear support vector machine is composed of a set of given support vectors z and a set of weights w. The computation for the output of a given SVM with N support vectors z₁, z₂, … , z_N and weights w₁, w₂, … , w_N is then given by:

$F(x) = sum_{i=1}^N w_i , left langle z_i,x right rangle + b$

Kernel Support Vector Machines

The original optimal hyperplane algorithm proposed by Vladimir Vapnik in 1963 was a linear classifier. However, in 1992, Bernhard Boser, Isabelle Guyon and Vapnik suggested a way to create non-linear classifiers by applying the kernel trick (originally proposed by Aizerman et al.) to maximum-margin hyperplanes. The resulting algorithm is formally similar, except that every dot product is replaced by a non-linear kernel function. This allows the algorithm to fit the maximum-margin hyperplane in a transformed feature space. The transformation may be non-linear and the transformed space high dimensional; thus though the classifier is a hyperplane in the high-dimensional feature space, it may be non-linear in the original input space.

Using kernels, the original formulation for the SVM given SVM with support vectors z₁, z₂, … , z_N and weights w₁, _w2, … , w_N is now given by:

$F(x) = sum_{i=1}^N w_i , k(z_i,x) + b$

It is also very straightforward to see that, using a linear kernel of the form K(z,x) = <z,x> = z^Tx, we recover the original formulation for the linear SVM.

The Kernel trick

The Kernel trick is a very interesting and powerful tool. It is powerful because it provides a bridge from linearity to non-linearity to any algorithm that solely depends on the dot product between two vectors. It comes from the fact that, if we first map our input data into a higher-dimensional space, a linear algorithm operating in this space will behave non-linearly in the original input space.

Now, the Kernel trick is really interesting because that mapping does not need to be ever computed. If our algorithm can be expressed only in terms of a inner product between two vectors, all we need is replace this inner product with the inner product from some other suitable space. That is where resides the “trick”: wherever a dot product is used, it is replaced with a Kernel function. The kernel function denotes an inner product in feature space and is usually denoted as:

K(x,y) = <φ(x),φ(y)>

Using the Kernel function, the algorithm can then be carried into a higher-dimension space without explicitly mapping the input points into this space. This is highly desirable, as sometimes our higher-dimensional feature space could even be infinite-dimensional and thus infeasible to compute.

Standard Kernels

Some common Kernel functions include the linear kernel, the polynomial kernel and the Gaussian kernel. Below is a simple list with their most interesting characteristics.

Linear Kernel	The Linear kernel is the simplest kernel function. It is given by the common inner product <x,y> plus an optional constant c. Kernel algorithms using a linear kernel are often equivalent to their non-kernel counterparts, i.e. KPCA with linear kernel is equivalent to standard PCA.	$k(x, y) = x^T y + c$
Polynomial Kernel	The Polynomial kernel is a non-stationary kernel. It is well suited for problems where all data is normalized.	$k(x, y) = (alpha x^T y + c)^d$
Gaussian Kernel	The Gaussian kernel is by far one of the most versatile Kernels. It is a radial basis function kernel, and is the preferred Kernel when we don’t know much about the data we are trying to model.	$k(x, y) = expleft(-frac{ lVert x-y rVert ^2}{2sigma^2}right)$

For more Kernel functions, check Kernel functions for Machine Learning Applications. The accompanying source code includes definitions for over 20 distinct Kernel functions, many of them detailed in the aforementioned post.

Learning Algorithms

Sequential Minimal Optimization

Previous SVM learning algorithms involved the use of quadratic programming solvers. Some of them used chunking to split the problem in smaller parts which could be solved more efficiently. Platt’s Sequential Minimal Optimization (SMO) algorithm puts chunking to the extreme by breaking the problem down into 2-dimensional sub-problems that can be solved analytically, eliminating the need for a numerical optimization algorithm.

The algorithm makes use of Lagrange multipliers to compute the optimization problem. Platt’s algorithm is composed of three main procedures or parts:

run, which iterates over all points until convergence to a tolerance threshold;
examineExample, which finds two points to jointly optimize;
takeStep, which solves the 2-dimensional optimization problem analytically.

The algorithm is also governed by three extra parameters besides the Kernel function and the data points.

The parameter Ccontrols the trade off between allowing some training errors and forcing rigid margins. Increasing the value of C increases the cost of misclassifications but may result in models that do not generalize well to points outside the training set.
The parameter ε controls the width of the ε-insensitive zone, used to fit the training data. The value of ε can affect the number of support vectors used to construct the regression function. The bigger ε, the fewer support vectors are selected and the solution becomes more sparse. On the other hand, increasing the ε-value by too much will result in less accurate models.
The parameter T is the convergence tolerance. It is the criterion for completing the training process.

After the algorithm ends, a new Support Vector Machine can be created using only the points whose Lagrange multipliers are higher than zero. The expected outputs y_i can be individually multiplied by their corresponding Lagrange multipliers a_i to form a single weight vector w.

$F(x) = sum_{i=0}^N { alpha_i y , k(z_i,x) } + b = sum_{i=0}^N { w_i , k(z_i,x) } + b$

Sequential Minimal Optimization for Regression

A version of SVM for regression was proposed in 1996 by Vladimir Vapnik, Harris Drucker, Chris Burges, Linda Kaufman and Alex Smola. The method was called support vector regression and, as is the case with the original Support Vector Machine formulation, depends only on a subset of the training data, because the cost function for building the model ignores any training data close to the model prediction that is within a tolerance threshold ε.

Platt’s algorithm has also been modified for regression. Albeit still maintaining much of its original structure, the difference lies in the fact that the modified algorithm uses two Lagrange multipliers â_i and a_i for each input point i. After the algorithm ends, a new Support Vector Machine can be created using only points whose both Lagrange multipliers are higher than zero. The multipliers â_i and a_i are then subtracted to form a single weight vector w.

$F(x) = sum_{i=0}^N { (hat{alpha_i} - alpha_i) , k(z_i,x) } + b = sum_{i=0}^N { w_i , k(z_i,x) } + b$

The algorithm is also governed by the same three parameters presented above. The parameter ε, however, receives a special meaning. It governs the size of the ε-insensitive tube over the regression line. The algorithm has been further developed and adapted by Alex J. Smola, Bernhard Schoelkopf and further optimizations were introduced by Shevade et al and Flake et al.

Source Code

Here is the class diagram for the Support Vector Machine module. We can see it is very simple in terms of standard class organization.

Class diagram for the (Kernel) Support Vector Machines module.

Support Vector Machine

Below is the class definition for the Linear Support Vector Machine. It is pretty much self explanatory.

/// <summary>
///   Sparse Linear Support Vector Machine (SVM)
/// </summary>
public class SupportVectorMachine
{
    private int inputCount;
    private double[][] supportVectors;
    private double[] weights;
    private double threshold;

    /// <summary>
    /// Creates a new Support Vector Machine
    /// </summary>
    public SupportVectorMachine(int inputs)
    {
        this.inputCount = inputs;
    }

    /// <summary>
    /// Gets the number of inputs accepted by this SVM.
    /// </summary>
    public int Inputs
    {
        get { return inputCount; }
    }

    /// <summary>
    /// Gets or sets the collection of support vectors used by this machine.
    /// </summary>
    public double[][] SupportVectors
    {
        get { return supportVectors; }
        set { supportVectors = value; }
    }

    /// <summary>
    /// Gets or sets the collection of weights used by this machine.
    /// </summary>
    public double[] Weights
    {
        get { return weights; }
        set { weights = value; }
    }

    /// <summary>
    /// Gets or sets the threshold (bias) term for this machine.
    /// </summary>
    public double Threshold
    {
        get { return threshold; }
        set { threshold = value; }
    }

    /// <summary>
    /// Computes the given input to produce the corresponding output.
    /// </summary>
    /// <param name="input">An input vector.</param>
    /// <returns>The ouput for the given input.</returns>
    public virtual double Compute(double[] input)
    {
        double s = threshold;
        for (int i = 0; i < supportVectors.Length; i++)
        {
            double p = 0;
            for (int j = 0; j < input.Length; j++)
                p += supportVectors[i][j] * input[j];
    
            s += weights[i] * p;
        }

        return s;
    }

    /// <summary>
    /// Computes the given inputs to produce the corresponding outputs.
    /// </summary>
    public double[] Compute(double[][] inputs)
    {
        double[] outputs = new double[inputs.Length];

        for (int i = 0; i < inputs.Length; i++)
            outputs[i] = Compute(inputs[i]);

        return outputs;
    }
}

/// <summary>

/// Sparse Linear Support Vector Machine (SVM)

/// </summary>

public class SupportVectorMachine

{

private int inputCount;

private double[][] supportVectors;

private double[] weights;

private double threshold;

/// <summary>

/// Creates a new Support Vector Machine

/// </summary>

public SupportVectorMachine(int inputs)

{

this.inputCount = inputs;

}

/// <summary>

/// Gets the number of inputs accepted by this SVM.

/// </summary>

public int Inputs

{

get { return inputCount; }

}

/// <summary>

/// Gets or sets the collection of support vectors used by this machine.

/// </summary>

public double[][] SupportVectors

{

get { return supportVectors; }

set { supportVectors = value; }

}

/// <summary>

/// Gets or sets the collection of weights used by this machine.

/// </summary>

public double[] Weights

{

get { return weights; }

set { weights = value; }

}

/// <summary>

/// Gets or sets the threshold (bias) term for this machine.

/// </summary>

public double Threshold

{

get { return threshold; }

set { threshold = value; }

}

/// <summary>

/// Computes the given input to produce the corresponding output.

/// </summary>

/// <param name="input">An input vector.</param>

/// <returns>The ouput for the given input.</returns>

public virtual double Compute(double[] input)

{

double s = threshold;

for (int i = 0; i < supportVectors.Length; i++)

{

double p = 0;

for (int j = 0; j < input.Length; j++)

p += supportVectors[i][j] * input[j];

s += weights[i] * p;

}

return s;

}

/// <summary>

/// Computes the given inputs to produce the corresponding outputs.

/// </summary>

public double[] Compute(double[][] inputs)

{

double[] outputs = new double[inputs.Length];

for (int i = 0; i < inputs.Length; i++)

outputs[i] = Compute(inputs[i]);

return outputs;

}

Kernel Support Vector Machine

Here is the class definition for the Kernel Support Vector Machine. It inherits from Support Vector Machine and extends it with a Kernel property. The Compute method is also overridden to include the chosen Kernel in the model computation.

/// <summary>
///  Sparse Kernel Support Vector Machine (kSVM)
/// </summary>
/// <remarks>
/// <para>
///   The original optimal hyperplane algorithm (SVM) proposed by Vladimir Vapnik in 1963 was a
///   linear classifier. However, in 1992, Bernhard Boser, Isabelle Guyon and Vapnik suggested
///   a way to create non-linear classifiers by applying the kernel trick (originally proposed
///   by Aizerman et al.) to maximum-margin hyperplanes. The resulting algorithm is formally
///   similar, except that every dot product is replaced by a non-linear kernel function.</para>
/// <para>
///   This allows the algorithm to fit the maximum-margin hyperplane in a transformed feature space.
///   The transformation may be non-linear and the transformed space high dimensional; thus though
///   the classifier is a hyperplane in the high-dimensional feature space, it may be non-linear in
///   the original input space.</para> 
/// <para>
///   References:
///   <list type="bullet">
///     <item><description><a href="http://en.wikipedia.org/wiki/Support_vector_machine">
///       http://en.wikipedia.org/wiki/Support_vector_machine</a></description></item>
///     <item><description><a href="http://www.kernel-machines.org/">
///       http://www.kernel-machines.org/</a></description></item>
///   </list></para>  
/// </remarks>
/// 
/// <example>
///   <code>
///   // Example XOR problem
///   double[][] inputs =
///   {
///       new double[] { 0, 0 }, // 0 xor 0: 1 (label +1)
///       new double[] { 0, 1 }, // 0 xor 1: 0 (label -1)
///       new double[] { 1, 0 }, // 1 xor 0: 0 (label -1)
///       new double[] { 1, 1 }  // 1 xor 1: 1 (label +1)
///   };
///   
///   // Dichotomy SVM outputs should be given as [-1;+1]
///   int[] labels =
///   {
///       // 1,  0,  0, 1
///          1, -1, -1, 1
///   };
///   
///   // Create a Kernel Support Vector Machine for the given inputs
///   KernelSupportVectorMachine machine = new KernelSupportVectorMachine(new Gaussian(0.1), inputs[0].Length);
///   
///   // Instantiate a new learning algorithm for SVMs
///   SequentialMinimalOptimization smo = new SequentialMinimalOptimization(machine, inputs, labels);
///   
///   // Set up the learning algorithm
///   smo.Complexity = 1.0;
///   
///   // Run the learning algorithm
///   double error = smo.Run();
///   
///   // Compute the decision output for one of the input vectors
///   int decision = System.Math.Sign(svm.Compute(inputs[0]));
///   </code>
/// </example>
/// 
[Serializable]
public class KernelSupportVectorMachine : SupportVectorMachine
{

    private IKernel kernel;


    /// <summary>
    ///   Creates a new Kernel Support Vector Machine.
    /// </summary>
    /// 
    /// <param name="kernel">The chosen kernel for the machine.</param>
    /// <param name="inputs">The number of inputs for the machine.</param>
    /// 
    /// <remarks>
    ///   If the number of inputs is zero, this means the machine
    ///   accepts a indefinite number of inputs. This is often the
    ///   case for kernel vector machines using a sequence kernel.
    /// </remarks>
    /// 
    public KernelSupportVectorMachine(IKernel kernel, int inputs) : base(inputs)
    {
        if (kernel == null) throw new ArgumentNullException("kernel");

        this.kernel = kernel;
    }

    /// <summary>
    ///   Gets or sets the kernel used by this machine.
    /// </summary>
    /// 
    public IKernel Kernel
    {
        get { return kernel; }
        set { kernel = value; }
    }

    /// <summary>
    ///   Computes the given input to produce the corresponding output.
    /// </summary>
    /// 
    /// <remarks>
    ///   For a binary decision problem, the decision for the negative
    ///   or positive class is typically computed by taking the sign of
    ///   the machine's output.
    /// </remarks>
    /// 
    /// <param name="inputs">An input vector.</param>
    /// <returns>The output for the given input.</returns>
    /// 
    public override double Compute(double[] inputs)
    {
        double s = Threshold;

        for (int i = 0; i < SupportVectors.Length; i++)
            s += Weights[i] * kernel.Function(SupportVectors[i], inputs);

        return s;
    }
}

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

/// <summary>

/// Sparse Kernel Support Vector Machine (kSVM)

/// </summary>

/// <remarks>

/// <para>

/// The original optimal hyperplane algorithm (SVM) proposed by Vladimir Vapnik in 1963 was a

/// linear classifier. However, in 1992, Bernhard Boser, Isabelle Guyon and Vapnik suggested

/// a way to create non-linear classifiers by applying the kernel trick (originally proposed

/// by Aizerman et al.) to maximum-margin hyperplanes. The resulting algorithm is formally

/// similar, except that every dot product is replaced by a non-linear kernel function.</para>

/// <para>

/// This allows the algorithm to fit the maximum-margin hyperplane in a transformed feature space.

/// The transformation may be non-linear and the transformed space high dimensional; thus though

/// the classifier is a hyperplane in the high-dimensional feature space, it may be non-linear in

/// the original input space.</para>

/// <para>

/// References:

/// <list type="bullet">

/// <item><description><a href="http://en.wikipedia.org/wiki/Support_vector_machine">

/// http://en.wikipedia.org/wiki/Support_vector_machine</a></description></item>

/// <item><description><a href="http://www.kernel-machines.org/">

/// http://www.kernel-machines.org/</a></description></item>

/// </list></para>

/// </remarks>

///

/// <example>

/// <code>

/// // Example XOR problem

/// double[][] inputs =

/// {

/// new double[] { 0, 0 }, // 0 xor 0: 1 (label +1)

/// new double[] { 0, 1 }, // 0 xor 1: 0 (label -1)

/// new double[] { 1, 0 }, // 1 xor 0: 0 (label -1)

/// new double[] { 1, 1 } // 1 xor 1: 1 (label +1)

/// };

///

/// // Dichotomy SVM outputs should be given as [-1;+1]

/// int[] labels =

/// {

/// // 1, 0, 0, 1

/// 1, -1, -1, 1

/// };

///

/// // Create a Kernel Support Vector Machine for the given inputs

/// KernelSupportVectorMachine machine = new KernelSupportVectorMachine(new Gaussian(0.1), inputs[0].Length);

///

/// // Instantiate a new learning algorithm for SVMs

/// SequentialMinimalOptimization smo = new SequentialMinimalOptimization(machine, inputs, labels);

///

/// // Set up the learning algorithm

/// smo.Complexity = 1.0;

///

/// // Run the learning algorithm

/// double error = smo.Run();

///

/// // Compute the decision output for one of the input vectors

/// int decision = System.Math.Sign(svm.Compute(inputs[0]));

/// </code>

/// </example>

///

[Serializable]

public class KernelSupportVectorMachine : SupportVectorMachine

{

private IKernel kernel;

/// <summary>

/// Creates a new Kernel Support Vector Machine.

/// </summary>

///

/// <param name="kernel">The chosen kernel for the machine.</param>

/// <param name="inputs">The number of inputs for the machine.</param>

///

/// <remarks>

/// If the number of inputs is zero, this means the machine

/// accepts a indefinite number of inputs. This is often the

/// case for kernel vector machines using a sequence kernel.

/// </remarks>

///

public KernelSupportVectorMachine(IKernel kernel, int inputs) : base(inputs)

{

if (kernel == null) throw new ArgumentNullException("kernel");

this.kernel = kernel;

}

/// <summary>

/// Gets or sets the kernel used by this machine.

/// </summary>

///

public IKernel Kernel

{

get { return kernel; }

set { kernel = value; }

}

/// <summary>

/// Computes the given input to produce the corresponding output.

/// </summary>

///

/// <remarks>

/// For a binary decision problem, the decision for the negative

/// or positive class is typically computed by taking the sign of

/// the machine's output.

/// </remarks>

///

/// <param name="inputs">An input vector.</param>

/// <returns>The output for the given input.</returns>

///

public override double Compute(double[] inputs)

{

double s = Threshold;

for (int i = 0; i < SupportVectors.Length; i++)

s += Weights[i] * kernel.Function(SupportVectors[i], inputs);

return s;

}

Sequential Minimal Optimization

Here is the code for the Sequential Minimal Optimization (SMO) algorithm.

/// <summary>
///   Sequential Minimal Optimization (SMO) Algorithm
/// </summary>
/// 
/// <remarks>
/// <para>
///   The SMO algorithm is an algorithm for solving large quadratic programming (QP)
///   optimization problems, widely used for the training of support vector machines.
///   First developed by John C. Platt in 1998, SMO breaks up large QP problems into
///   a series of smallest possible QP problems, which are then solved analytically.</para>
/// <para>
///   This class follows the original algorithm by Platt as strictly as possible.</para>
///  
/// <para>
///   References:
///   <list type="bullet">
///     <item><description>
///       <a href="http://en.wikipedia.org/wiki/Sequential_Minimal_Optimization">
///       Wikipedia, The Free Encyclopedia. Sequential Minimal Optimization. Available on:
///       http://en.wikipedia.org/wiki/Sequential_Minimal_Optimization </a></description></item>
///     <item><description>
///       <a href="http://research.microsoft.com/en-us/um/people/jplatt/smoTR.pdf">
///       John C. Platt, Sequential Minimal Optimization: A Fast Algorithm for Training Support
///       Vector Machines. 1998. Available on: http://research.microsoft.com/en-us/um/people/jplatt/smoTR.pdf </a></description></item>
///     <item><description>
///       <a href="http://www.idiom.com/~zilla/Work/Notes/svmtutorial.pdf">
///       J. P. Lewis. A Short SVM (Support Vector Machine) Tutorial. Available on:
///       http://www.idiom.com/~zilla/Work/Notes/svmtutorial.pdf </a></description></item>
///     </list></para>  
/// </remarks>
/// 
/// <example>
///   <code>
///   // Example XOR problem
///   double[][] inputs =
///   {
///       new double[] { 0, 0 }, // 0 xor 0: 1 (label +1)
///       new double[] { 0, 1 }, // 0 xor 1: 0 (label -1)
///       new double[] { 1, 0 }, // 1 xor 0: 0 (label -1)
///       new double[] { 1, 1 }  // 1 xor 1: 1 (label +1)
///   };
///    
///   // Dichotomy SVM outputs should be given as [-1;+1]
///   int[] labels =
///   {
///          1, -1, -1, 1
///   };
///  
///   // Create a Kernel Support Vector Machine for the given inputs
///   KernelSupportVectorMachine machine = new KernelSupportVectorMachine(new Gaussian(0.1), inputs[0].Length);
/// 
///   // Instantiate a new learning algorithm for SVMs
///   SequentialMinimalOptimization smo = new SequentialMinimalOptimization(machine, inputs, labels);
/// 
///   // Set up the learning algorithm
///   smo.Complexity = 1.0;
/// 
///   // Run the learning algorithm
///   double error = smo.Run();
///   
///   // Compute the decision output for one of the input vectors
///   int decision = System.Math.Sign(svm.Compute(inputs[0]));
///  </code>
/// </example>
/// 
public class SequentialMinimalOptimization : ISupportVectorMachineLearning
{
    private static Random random = new Random();

    // Training data
    private double[][] inputs;
    private int[] outputs;

    // Learning algorithm parameters
    private double c = 1.0;
    private double tolerance = 1e-3;
    private double epsilon = 1e-3;
    private bool useComplexityHeuristic;

    // Support Vector Machine parameters
    private SupportVectorMachine machine;
    private IKernel kernel;
    private double[] alpha;
    private double bias;


    // Error cache to speed up computations
    private double[] errors;


    /// <summary>
    ///   Initializes a new instance of a Sequential Minimal Optimization (SMO) algorithm.
    /// </summary>
    /// 
    /// <param name="machine">A Support Vector Machine.</param>
    /// <param name="inputs">The input data points as row vectors.</param>
    /// <param name="outputs">The classification label for each data point in the range [-1;+1].</param>
    /// 
    public SequentialMinimalOptimization(SupportVectorMachine machine,
        double[][] inputs, int[] outputs)
    {

        // Initial argument checking
        if (machine == null)
            throw new ArgumentNullException("machine");

        if (inputs == null)
            throw new ArgumentNullException("inputs");

        if (outputs == null)
            throw new ArgumentNullException("outputs");

        if (inputs.Length != outputs.Length)
            throw new ArgumentException("The number of inputs and outputs does not match.", "outputs");

        for (int i = 0; i < outputs.Length; i++)
        {
            if (outputs[i] != 1 && outputs[i] != -1)
                throw new ArgumentOutOfRangeException("outputs", "One of the labels in the output vector is neither +1 or -1.");
        }

        if (machine.Inputs > 0)
        {
            // This machine has a fixed input vector size
            for (int i = 0; i < inputs.Length; i++)
                if (inputs[i].Length != machine.Inputs)
                    throw new ArgumentException("The size of the input vectors does not match the expected number of inputs of the machine");
        }


        // Machine
        this.machine = machine;

        // Kernel (if applicable)
        KernelSupportVectorMachine ksvm = machine as KernelSupportVectorMachine;
        this.kernel = (ksvm != null) ? ksvm.Kernel : new Linear();


        // Learning data
        this.inputs = inputs;
        this.outputs = outputs;

    }


    //---------------------------------------------


    #region Properties
    /// <summary>
    ///   Complexity (cost) parameter C. Increasing the value of C forces the creation
    ///   of a more accurate model that may not generalize well. Default value is the
    ///   number of examples divided by the trace of the kernel matrix.
    /// </summary>
    /// <remarks>
    ///   The cost parameter C controls the trade off between allowing training
    ///   errors and forcing rigid margins. It creates a soft margin that permits
    ///   some misclassifications. Increasing the value of C increases the cost of
    ///   misclassifying points and forces the creation of a more accurate model
    ///   that may not generalize well.
    /// </remarks>
    public double Complexity
    {
        get { return this.c; }
        set { this.c = value; }
    }


    /// <summary>
    ///   Gets or sets a value indicating whether the Complexity parameter C
    ///   should be computed automatically by employing an heuristic rule.
    /// </summary>
    /// <value>
    ///     <c>true</c> if complexity should be computed automatically; otherwise, <c>false</c>.
    /// </value>
    public bool UseComplexityHeuristic
    {
        get { return useComplexityHeuristic; }
        set { useComplexityHeuristic = value; }
    }

    /// <summary>
    ///   Insensitivity zone ε. Increasing the value of ε can result in fewer support
    ///   vectors in the created model. Default value is 1e-3.
    /// </summary>
    /// <remarks>
    ///   Parameter ε controls the width of the ε-insensitive zone, used to fit the training
    ///   data. The value of ε can affect the number of support vectors used to construct the
    ///   regression function. The bigger ε, the fewer support vectors are selected. On the
    ///   other hand, bigger ε-values results in more flat estimates.
    /// </remarks>
    public double Epsilon
    {
        get { return epsilon; }
        set { epsilon = value; }
    }

    /// <summary>
    ///   Convergence tolerance. Default value is 1e-3.
    /// </summary>
    /// <remarks>
    ///   The criterion for completing the model training process. The default is 0.001.
    /// </remarks>
    public double Tolerance
    {
        get { return this.tolerance; }
        set { this.tolerance = value; }
    }
    #endregion


    //---------------------------------------------


    /// <summary>
    ///   Runs the SMO algorithm.
    /// </summary>
    /// 
    /// <param name="computeError">
    ///   True to compute error after the training
    ///   process completes, false otherwise. Default is true.
    /// </param>
    /// 
    /// <returns>
    ///   The misclassification error rate of
    ///   the resulting support vector machine.
    /// </returns>
    /// 
    public double Run(bool computeError)
    {

        // The SMO algorithm chooses to solve the smallest possible optimization problem
        // at every step. At every step, SMO chooses two Lagrange multipliers to jointly
        // optimize, finds the optimal values for these multipliers, and updates the SVM
        // to reflect the new optimal values
        //
        // Reference: http://research.microsoft.com/en-us/um/people/jplatt/smoTR.pdf


        // Initialize variables
        int N = inputs.Length;

        if (useComplexityHeuristic)
            c = computeComplexity();

        // Lagrange multipliers
        this.alpha = new double[N];

        // Error cache
        this.errors = new double[N];

        // Algorithm:
        int numChanged = 0;
        int examineAll = 1;

        while (numChanged > 0 || examineAll > 0)
        {
            numChanged = 0;
            if (examineAll > 0)
            {
                // loop I over all training examples
                for (int i = 0; i < N; i++)
                    numChanged += examineExample(i);
            }
            else
            {
                // loop I over examples where alpha is not 0 and not C
                for (int i = 0; i < N; i++)
                    if (alpha[i] != 0 && alpha[i] != c)
                        numChanged += examineExample(i);
            }

            if (examineAll == 1)
                examineAll = 0;
            else if (numChanged == 0)
                examineAll = 1;
        }


        // Store Support Vectors in the SV Machine. Only vectors which have lagrange multipliers
        // greater than zero will be stored as only those are actually required during evaluation.
        List<int> indices = new List<int>();
        for (int i = 0; i < N; i++)
        {
            // Only store vectors with multipliers > 0
            if (alpha[i] > 0) indices.Add(i);
        }

        int vectors = indices.Count;
        machine.SupportVectors = new double[vectors][];
        machine.Weights = new double[vectors];
        for (int i = 0; i < vectors; i++)
        {
            int j = indices[i];
            machine.SupportVectors[i] = inputs[j];
            machine.Weights[i] = alpha[j] * outputs[j];
        }
        machine.Threshold = -bias;


        // Compute error if required.
        return (computeError) ? ComputeError(inputs, outputs) : 0.0;
    }

    /// <summary>
    ///   Runs the SMO algorithm.
    /// </summary>
    /// 
    /// <returns>
    ///   The misclassification error rate of
    ///   the resulting support vector machine.
    /// </returns>
    /// 
    public double Run()
    {
        return Run(true);
    }

    /// <summary>
    ///   Computes the error rate for a given set of input and outputs.
    /// </summary>
    /// 
    public double ComputeError(double[][] inputs, int[] expectedOutputs)
    {
        // Compute errors
        int count = 0;
        for (int i = 0; i < inputs.Length; i++)
        {
            if (Math.Sign(compute(inputs[i])) != Math.Sign(expectedOutputs[i]))
                count++;
        }

        // Return misclassification error ratio
        return (double)count / inputs.Length;
    }

    //---------------------------------------------


    /// <summary>
    ///  Chooses which multipliers to optimize using heuristics.
    /// </summary>
    /// 
    private int examineExample(int i2)
    {
        double[] p2 = inputs[i2]; // Input point at index i2
        double y2 = outputs[i2];  // Classification label for p2
        double alph2 = alpha[i2];    // Lagrange multiplier for p2

        // SVM output on p2 - y2. Check if it has already been computed
        double e2 = (alph2 > 0 && alph2 < c) ? errors[i2] : compute(p2) - y2;

        double r2 = y2 * e2;


        // Heuristic 01 (for the first multiplier choice):
        //  - Testing for KKT conditions within the tolerance margin
        if (!(r2 < -tolerance && alph2 < c) && !(r2 > tolerance && alph2 > 0))
            return 0;


        // Heuristic 02 (for the second multiplier choice):
        //  - Once a first Lagrange multiplier is chosen, SMO chooses the second Lagrange multiplier to
        //    maximize the size of the step taken during joint optimization. Now, evaluating the kernel
        //    function is time consuming, so SMO approximates the step size by the absolute value of the
        //    absolute error difference.
        int i1 = -1; double max = 0;
        for (int i = 0; i < inputs.Length; i++)
        {
            if (alpha[i] > 0 && alpha[i] < c)
            {
                double error1 = errors[i];
                double aux = System.Math.Abs(e2 - error1);

                if (aux > max)
                {
                    max = aux;
                    i1 = i;
                }
            }
        }

        if (i1 >= 0 && takeStep(i1, i2)) return 1;


        // Heuristic 03:
        //  - Under unusual circumstances, SMO cannot make positive progress using the second
        //    choice heuristic above. If it is the case, then SMO starts iterating through the
        //    non-bound examples, searching for an second example that can make positive progress.

        int start = random.Next(inputs.Length);
        for (i1 = start; i1 < inputs.Length; i1++)
        {
            if (alpha[i1] > 0 && alpha[i1] < c)
                if (takeStep(i1, i2)) return 1;
        }
        for (i1 = 0; i1 < start; i1++)
        {
            if (alpha[i1] > 0 && alpha[i1] < c)
                if (takeStep(i1, i2)) return 1;
        }


        // Heuristic 04:
        //  - If none of the non-bound examples make positive progress, then SMO starts iterating
        //    through the entire training set until an example is found that makes positive progress.
        //    Both the iteration through the non-bound examples and the iteration through the entire
        //    training set are started at random locations, in order not to bias SMO towards the
        //    examples at the beginning of the training set. 

        start = random.Next(inputs.Length);
        for (i1 = start; i1 < inputs.Length; i1++)
        {
            if (takeStep(i1, i2)) return 1;
        }
        for (i1 = 0; i1 < start; i1++)
        {
            if (takeStep(i1, i2)) return 1;
        }


        // In extremely degenerate circumstances, none of the examples will make an adequate second
        // example. When this happens, the first example is skipped and SMO continues with another
        // chosen first example.
        return 0;
    }

    /// <summary>
    ///   Analytically solves the optimization problem for two Lagrange multipliers.
    /// </summary>
    /// 
    private bool takeStep(int i1, int i2)
    {
        if (i1 == i2) return false;

        double[] p1 = inputs[i1]; // Input point at index i1
        double alph1 = alpha[i1];    // Lagrange multiplier for p1
        double y1 = outputs[i1];  // Classification label for p1

        // SVM output on p1 - y1. Check if it has already been computed
        double e1 = (alph1 > 0 && alph1 < c) ? errors[i1] : compute(p1) - y1;

        double[] p2 = inputs[i2]; // Input point at index i2
        double alph2 = alpha[i2];    // Lagrange multiplier for p2
        double y2 = outputs[i2];  // Classification label for p2

        // SVM output on p2 - y2. Check if it has already been computed
        double e2 = (alph2 > 0 && alph2 < c) ? errors[i2] : compute(p2) - y2;


        double s = y1 * y2;


        // Compute L and H according to equations (13) and (14) (Platt, 1998)
        double L, H;
        if (y1 != y2)
        {
            // If the target y1 does not equal the target           (13)
            // y2, then the following bounds apply to a2:
            L = Math.Max(0, alph2 - alph1);
            H = Math.Min(c, c + alph2 - alph1);
        }
        else
        {
            // If the target y1 does equal the target               (14)
            // y2, then the following bounds apply to a2:
            L = Math.Max(0,  alph2 + alph1 - c);
            H = Math.Min(c,  alph2 + alph1);
        }

        if (L == H) return false;

        double k11, k22, k12, eta;
        k11 = kernel.Function(p1, p1);
        k12 = kernel.Function(p1, p2);
        k22 = kernel.Function(p2, p2);
        eta = k11 + k22 - 2.0 * k12;

        double a1, a2;

        if (eta > 0)
        {
            a2 = alph2 - y2 * (e2 - e1) / eta;

            if (a2 < L) a2 = L;
            else if (a2 > H) a2 = H;
        }
        else
        {
            // Compute objective function Lobj and Hobj at
            //  a2=L and a2=H respectively, using (eq. 19)

            double L1 = alph1 + s * (alph2 - L);
            double H1 = alph1 + s * (alph2 - H);
            double f1 = y1 * (e1 + bias) - alph1 * k11 - s * alph2 * k12;
            double f2 = y2 * (e2 + bias) - alph2 * k22 - s * alph1 * k12;
            double Lobj = -0.5 * L1 * L1 * k11 - 0.5 * L * L * k22 - s * L * L1 * k12 - L1 * f1 - L * f2;
            double Hobj = -0.5 * H1 * H1 * k11 - 0.5 * H * H * k22 - s * H * H1 * k12 - H1 * f1 - H * f2;

            if (Lobj > Hobj + epsilon) a2 = L;
            else if (Lobj < Hobj - epsilon) a2 = H;
            else a2 = alph2;
        }

        if (Math.Abs(a2 - alph2) < epsilon * (a2 + alph2 + epsilon))
            return false;

        a1 = alph1 + s * (alph2 - a2);

        if (a1 < 0)
        {
            a2 += s * a1;
            a1 = 0;
        }
        else if (a1 > c)
        {
            double d = a1 - c;
            a2 += s * d;
            a1 = c;
        }


        // Update threshold (bias) to reflect change in Lagrange multipliers
        double b1 = 0, b2 = 0;
        double new_b = 0, delta_b;

        if (a1 > 0 && a1 < c)
        {
            // a1 is not at bounds
            new_b = e1 + y1 * (a1 - alph1) * k11 + y2 * (a2 - alph2) * k12 + bias;
        }
        else
        {
            if (a2 > 0 && a2 < c)
            {
                // a1 is at bounds but a2 is not.
                new_b = e2 + y1 * (a1 - alph1) * k12 + y2 * (a2 - alph2) * k22 + bias;
            }
            else
            {
                // Both new Lagrange multipliers are at bound. SMO algorithm
                // chooses the threshold to be halfway in between b1 and b2.
                b1 = e1 + y1 * (a1 - alph1) * k11 + y2 * (a2 - alph2) * k12 + bias;
                b2 = e2 + y1 * (a1 - alph1) * k12 + y2 * (a2 - alph2) * k22 + bias;
                new_b = (b1 + b2) / 2;
            }
        }

        delta_b = new_b - bias;
        bias = new_b;


        // Update error cache using new Lagrange multipliers
        double t1 = y1 * (a1 - alph1);
        double t2 = y2 * (a2 - alph2);

        for (int i = 0; i < inputs.Length; i++)
        {
            if (0 < alpha[i] && alpha[i] < c)
            {
                double[] point = inputs[i];
                errors[i] +=
                        t1 * kernel.Function(p1, point) +
                        t2 * kernel.Function(p2, point) -
                        delta_b;
            }
        }

        errors[i1] = 0f;
        errors[i2] = 0f;


        // Update lagrange multipliers
        alpha[i1] = a1;
        alpha[i2] = a2;


        return true;
    }

    /// <summary>
    ///   Computes the SVM output for a given point.
    /// </summary>
    /// 
    private double compute(double[] point)
    {
        double sum = -bias;
        for (int i = 0; i < inputs.Length; i++)
        {
            if (alpha[i] > 0)
                sum += alpha[i] * outputs[i] * kernel.Function(inputs[i], point);
        }

        return sum;
    }


    private double computeComplexity()
    {
        // Compute initial value for C as the number of examples
        // divided by the trace of the input sample kernel matrix.
        double sum = 0.0;
        for (int i = 0; i < inputs.Length; i++)
            sum += kernel.Function(inputs[i], inputs[i]);
        return inputs.Length / sum;
    }
}

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300

301

302

303

304

305

306

307

308

309

310

311

312

313

314

315

316

317

318

319

320

321

322

323

324

325

326

327

328

329

330

331

332

333

334

335

336

337

338

339

340

341

342

343

344

345

346

347

348

349

350

351

352

353

354

355

356

357

358

359

360

361

362

363

364

365

366

367

368

369

370

371

372

373

374

375

376

377

378

379

380

381

382

383

384

385

386

387

388

389

390

391

392

393

394

395

396

397

398

399

400

401

402

403

404

405

406

407

408

409

410

411

412

413

414

415

416

417

418

419

420

421

422

423

424

425

426

427

428

429

430

431

432

433

434

435

436

437

438

439

440

441

442

443

444

445

446

447

448

449

450

451

452

453

454

455

456

457

458

459

460

461

462

463

464

465

466

467

468

469

470

471

472

473

474

475

476

477

478

479

480

481

482

483

484

485

486

487

488

489

490

491

492

493

494

495

496

497

498

499

500

501

502

503

504

505

506

507

508

509

510

511

512

513

514

515

516

517

518

519

520

521

522

523

524

525

526

527

528

529

530

531

532

533

534

535

536

537

538

539

540

541

542

543

544

545

546

547

548

549

550

551

552

553

554

555

556

557

558

559

560

561

562

563

564

565

566

567

568

569

570

571

572

573

574

575

576

577

578

579

580

581

582

583

584

585

586

587

588

589

590

591

592

593

594

595

596

597

598

599

600

601

602

603

604

605

606

607

/// <summary>

/// Sequential Minimal Optimization (SMO) Algorithm

/// </summary>

///

/// <remarks>

/// <para>

/// The SMO algorithm is an algorithm for solving large quadratic programming (QP)

/// optimization problems, widely used for the training of support vector machines.

/// First developed by John C. Platt in 1998, SMO breaks up large QP problems into

/// a series of smallest possible QP problems, which are then solved analytically.</para>

/// <para>

/// This class follows the original algorithm by Platt as strictly as possible.</para>

///

/// <para>

/// References:

/// <list type="bullet">

/// <item><description>

/// <a href="http://en.wikipedia.org/wiki/Sequential_Minimal_Optimization">

/// Wikipedia, The Free Encyclopedia. Sequential Minimal Optimization. Available on:

/// http://en.wikipedia.org/wiki/Sequential_Minimal_Optimization </a></description></item>

/// <item><description>

/// <a href="http://research.microsoft.com/en-us/um/people/jplatt/smoTR.pdf">

/// John C. Platt, Sequential Minimal Optimization: A Fast Algorithm for Training Support

/// Vector Machines. 1998. Available on: http://research.microsoft.com/en-us/um/people/jplatt/smoTR.pdf </a></description></item>

/// <item><description>

/// <a href="http://www.idiom.com/~zilla/Work/Notes/svmtutorial.pdf">

/// J. P. Lewis. A Short SVM (Support Vector Machine) Tutorial. Available on:

/// http://www.idiom.com/~zilla/Work/Notes/svmtutorial.pdf </a></description></item>

/// </list></para>

/// </remarks>

///

/// <example>

/// <code>

/// // Example XOR problem

/// double[][] inputs =

/// {

/// new double[] { 0, 0 }, // 0 xor 0: 1 (label +1)

/// new double[] { 0, 1 }, // 0 xor 1: 0 (label -1)

/// new double[] { 1, 0 }, // 1 xor 0: 0 (label -1)

/// new double[] { 1, 1 } // 1 xor 1: 1 (label +1)

/// };

///

/// // Dichotomy SVM outputs should be given as [-1;+1]

/// int[] labels =

/// {

/// 1, -1, -1, 1

/// };

///

/// // Create a Kernel Support Vector Machine for the given inputs

/// KernelSupportVectorMachine machine = new KernelSupportVectorMachine(new Gaussian(0.1), inputs[0].Length);

///

/// // Instantiate a new learning algorithm for SVMs

/// SequentialMinimalOptimization smo = new SequentialMinimalOptimization(machine, inputs, labels);

///

/// // Set up the learning algorithm

/// smo.Complexity = 1.0;

///

/// // Run the learning algorithm

/// double error = smo.Run();

///

/// // Compute the decision output for one of the input vectors

/// int decision = System.Math.Sign(svm.Compute(inputs[0]));

/// </code>

/// </example>

///

public class SequentialMinimalOptimization : ISupportVectorMachineLearning

{

private static Random random = new Random();

// Training data

private double[][] inputs;

private int[] outputs;

// Learning algorithm parameters

private double c = 1.0;

private double tolerance = 1e-3;

private double epsilon = 1e-3;

private bool useComplexityHeuristic;

// Support Vector Machine parameters

private SupportVectorMachine machine;

private IKernel kernel;

private double[] alpha;

private double bias;

// Error cache to speed up computations

private double[] errors;

/// <summary>

/// Initializes a new instance of a Sequential Minimal Optimization (SMO) algorithm.

/// </summary>

///

/// <param name="machine">A Support Vector Machine.</param>

/// <param name="inputs">The input data points as row vectors.</param>

/// <param name="outputs">The classification label for each data point in the range [-1;+1].</param>

///

public SequentialMinimalOptimization(SupportVectorMachine machine,

double[][] inputs, int[] outputs)

{

// Initial argument checking

if (machine == null)

throw new ArgumentNullException("machine");

if (inputs == null)

throw new ArgumentNullException("inputs");

if (outputs == null)

throw new ArgumentNullException("outputs");

if (inputs.Length != outputs.Length)

throw new ArgumentException("The number of inputs and outputs does not match.", "outputs");

for (int i = 0; i < outputs.Length; i++)

{

if (outputs[i] != 1 && outputs[i] != -1)

throw new ArgumentOutOfRangeException("outputs", "One of the labels in the output vector is neither +1 or -1.");

}

if (machine.Inputs > 0)

{

// This machine has a fixed input vector size

for (int i = 0; i < inputs.Length; i++)

if (inputs[i].Length != machine.Inputs)

throw new ArgumentException("The size of the input vectors does not match the expected number of inputs of the machine");

}

// Machine

this.machine = machine;

// Kernel (if applicable)

KernelSupportVectorMachine ksvm = machine as KernelSupportVectorMachine;

this.kernel = (ksvm != null) ? ksvm.Kernel : new Linear();

// Learning data

this.inputs = inputs;

this.outputs = outputs;

}

//---------------------------------------------

#region Properties

/// <summary>

/// Complexity (cost) parameter C. Increasing the value of C forces the creation

/// of a more accurate model that may not generalize well. Default value is the

/// number of examples divided by the trace of the kernel matrix.

/// </summary>

/// <remarks>

/// The cost parameter C controls the trade off between allowing training

/// errors and forcing rigid margins. It creates a soft margin that permits

/// some misclassifications. Increasing the value of C increases the cost of

/// misclassifying points and forces the creation of a more accurate model

/// that may not generalize well.

/// </remarks>

public double Complexity

{

get { return this.c; }

set { this.c = value; }

}

/// <summary>

/// Gets or sets a value indicating whether the Complexity parameter C

/// should be computed automatically by employing an heuristic rule.

/// </summary>

/// <value>

/// <c>true</c> if complexity should be computed automatically; otherwise, <c>false</c>.

/// </value>

public bool UseComplexityHeuristic

{

get { return useComplexityHeuristic; }

set { useComplexityHeuristic = value; }

}

/// <summary>

/// Insensitivity zone ε. Increasing the value of ε can result in fewer support

/// vectors in the created model. Default value is 1e-3.

/// </summary>

/// <remarks>

/// Parameter ε controls the width of the ε-insensitive zone, used to fit the training

/// data. The value of ε can affect the number of support vectors used to construct the

/// regression function. The bigger ε, the fewer support vectors are selected. On the

/// other hand, bigger ε-values results in more flat estimates.

/// </remarks>

public double Epsilon

{

get { return epsilon; }

set { epsilon = value; }

}

/// <summary>

/// Convergence tolerance. Default value is 1e-3.

/// </summary>

/// <remarks>

/// The criterion for completing the model training process. The default is 0.001.

/// </remarks>

public double Tolerance

{

get { return this.tolerance; }

set { this.tolerance = value; }

}

#endregion

//---------------------------------------------

/// <summary>

/// Runs the SMO algorithm.

/// </summary>

///

/// <param name="computeError">

/// True to compute error after the training

/// process completes, false otherwise. Default is true.

/// </param>

///

/// <returns>

/// The misclassification error rate of

/// the resulting support vector machine.

/// </returns>

///

public double Run(bool computeError)

{

// The SMO algorithm chooses to solve the smallest possible optimization problem

// at every step. At every step, SMO chooses two Lagrange multipliers to jointly

// optimize, finds the optimal values for these multipliers, and updates the SVM

// to reflect the new optimal values

// Reference: http://research.microsoft.com/en-us/um/people/jplatt/smoTR.pdf

// Initialize variables

int N = inputs.Length;

if (useComplexityHeuristic)

c = computeComplexity();

// Lagrange multipliers

this.alpha = new double[N];

// Error cache

this.errors = new double[N];

// Algorithm:

int numChanged = 0;

int examineAll = 1;

while (numChanged > 0 || examineAll > 0)

{

numChanged = 0;

if (examineAll > 0)

{

// loop I over all training examples

for (int i = 0; i < N; i++)

numChanged += examineExample(i);

}

else

{

// loop I over examples where alpha is not 0 and not C

for (int i = 0; i < N; i++)

if (alpha[i] != 0 && alpha[i] != c)

numChanged += examineExample(i);

}

if (examineAll == 1)

examineAll = 0;

else if (numChanged == 0)

examineAll = 1;

}

// Store Support Vectors in the SV Machine. Only vectors which have lagrange multipliers

// greater than zero will be stored as only those are actually required during evaluation.

List<int> indices = new List<int>();

for (int i = 0; i < N; i++)

{

// Only store vectors with multipliers > 0

if (alpha[i] > 0) indices.Add(i);

}

int vectors = indices.Count;

machine.SupportVectors = new double[vectors][];

machine.Weights = new double[vectors];

for (int i = 0; i < vectors; i++)

{

int j = indices[i];

machine.SupportVectors[i] = inputs[j];

machine.Weights[i] = alpha[j] * outputs[j];

}

machine.Threshold = -bias;

// Compute error if required.

return (computeError) ? ComputeError(inputs, outputs) : 0.0;

}

/// <summary>

/// Runs the SMO algorithm.

/// </summary>

///

/// <returns>

/// The misclassification error rate of

/// the resulting support vector machine.

/// </returns>

///

public double Run()

{

return Run(true);

}

/// <summary>

/// Computes the error rate for a given set of input and outputs.

/// </summary>

///

public double ComputeError(double[][] inputs, int[] expectedOutputs)

{

// Compute errors

int count = 0;

for (int i = 0; i < inputs.Length; i++)

{

if (Math.Sign(compute(inputs[i])) != Math.Sign(expectedOutputs[i]))

count++;

}

// Return misclassification error ratio

return (double)count / inputs.Length;

}

//---------------------------------------------

/// <summary>

/// Chooses which multipliers to optimize using heuristics.

/// </summary>

///

private int examineExample(int i2)

{

double[] p2 = inputs[i2]; // Input point at index i2

double y2 = outputs[i2]; // Classification label for p2

double alph2 = alpha[i2]; // Lagrange multiplier for p2

// SVM output on p2 - y2. Check if it has already been computed

double e2 = (alph2 > 0 && alph2 < c) ? errors[i2] : compute(p2) - y2;

double r2 = y2 * e2;

// Heuristic 01 (for the first multiplier choice):

// - Testing for KKT conditions within the tolerance margin

if (!(r2 < -tolerance && alph2 < c) && !(r2 > tolerance && alph2 > 0))

return 0;

// Heuristic 02 (for the second multiplier choice):

// - Once a first Lagrange multiplier is chosen, SMO chooses the second Lagrange multiplier to

// maximize the size of the step taken during joint optimization. Now, evaluating the kernel

// function is time consuming, so SMO approximates the step size by the absolute value of the

// absolute error difference.

int i1 = -1; double max = 0;

for (int i = 0; i < inputs.Length; i++)

{

if (alpha[i] > 0 && alpha[i] < c)

{

double error1 = errors[i];

double aux = System.Math.Abs(e2 - error1);

if (aux > max)

{

max = aux;

i1 = i;

}

if (i1 >= 0 && takeStep(i1, i2)) return 1;

// Heuristic 03:

// - Under unusual circumstances, SMO cannot make positive progress using the second

// choice heuristic above. If it is the case, then SMO starts iterating through the

// non-bound examples, searching for an second example that can make positive progress.

int start = random.Next(inputs.Length);

for (i1 = start; i1 < inputs.Length; i1++)

{

if (alpha[i1] > 0 && alpha[i1] < c)

if (takeStep(i1, i2)) return 1;

}

for (i1 = 0; i1 < start; i1++)

{

if (alpha[i1] > 0 && alpha[i1] < c)

if (takeStep(i1, i2)) return 1;

}

// Heuristic 04:

// - If none of the non-bound examples make positive progress, then SMO starts iterating

// through the entire training set until an example is found that makes positive progress.

// Both the iteration through the non-bound examples and the iteration through the entire

// training set are started at random locations, in order not to bias SMO towards the

// examples at the beginning of the training set.

start = random.Next(inputs.Length);

for (i1 = start; i1 < inputs.Length; i1++)

{

if (takeStep(i1, i2)) return 1;

}

for (i1 = 0; i1 < start; i1++)

{

if (takeStep(i1, i2)) return 1;

}

// In extremely degenerate circumstances, none of the examples will make an adequate second

// example. When this happens, the first example is skipped and SMO continues with another

// chosen first example.

return 0;

}

/// <summary>

/// Analytically solves the optimization problem for two Lagrange multipliers.

/// </summary>

///

private bool takeStep(int i1, int i2)

{

if (i1 == i2) return false;

double[] p1 = inputs[i1]; // Input point at index i1

double alph1 = alpha[i1]; // Lagrange multiplier for p1

double y1 = outputs[i1]; // Classification label for p1

// SVM output on p1 - y1. Check if it has already been computed

double e1 = (alph1 > 0 && alph1 < c) ? errors[i1] : compute(p1) - y1;

double[] p2 = inputs[i2]; // Input point at index i2

double alph2 = alpha[i2]; // Lagrange multiplier for p2

double y2 = outputs[i2]; // Classification label for p2

// SVM output on p2 - y2. Check if it has already been computed

double e2 = (alph2 > 0 && alph2 < c) ? errors[i2] : compute(p2) - y2;

double s = y1 * y2;

// Compute L and H according to equations (13) and (14) (Platt, 1998)

double L, H;

if (y1 != y2)

{

// If the target y1 does not equal the target (13)

// y2, then the following bounds apply to a2:

L = Math.Max(0, alph2 - alph1);

H = Math.Min(c, c + alph2 - alph1);

}

else

{

// If the target y1 does equal the target (14)

// y2, then the following bounds apply to a2:

L = Math.Max(0, alph2 + alph1 - c);

H = Math.Min(c, alph2 + alph1);

}

if (L == H) return false;

double k11, k22, k12, eta;

k11 = kernel.Function(p1, p1);

k12 = kernel.Function(p1, p2);

k22 = kernel.Function(p2, p2);

eta = k11 + k22 - 2.0 * k12;

double a1, a2;

if (eta > 0)

{

a2 = alph2 - y2 * (e2 - e1) / eta;

if (a2 < L) a2 = L;

else if (a2 > H) a2 = H;

}

else

{

// Compute objective function Lobj and Hobj at

// a2=L and a2=H respectively, using (eq. 19)

double L1 = alph1 + s * (alph2 - L);

double H1 = alph1 + s * (alph2 - H);

double f1 = y1 * (e1 + bias) - alph1 * k11 - s * alph2 * k12;

double f2 = y2 * (e2 + bias) - alph2 * k22 - s * alph1 * k12;

double Lobj = -0.5 * L1 * L1 * k11 - 0.5 * L * L * k22 - s * L * L1 * k12 - L1 * f1 - L * f2;

double Hobj = -0.5 * H1 * H1 * k11 - 0.5 * H * H * k22 - s * H * H1 * k12 - H1 * f1 - H * f2;

if (Lobj > Hobj + epsilon) a2 = L;

else if (Lobj < Hobj - epsilon) a2 = H;

else a2 = alph2;

}

if (Math.Abs(a2 - alph2) < epsilon * (a2 + alph2 + epsilon))

return false;

a1 = alph1 + s * (alph2 - a2);

if (a1 < 0)

{

a2 += s * a1;

a1 = 0;

}

else if (a1 > c)

{

double d = a1 - c;

a2 += s * d;

a1 = c;

}

// Update threshold (bias) to reflect change in Lagrange multipliers

double b1 = 0, b2 = 0;

double new_b = 0, delta_b;

if (a1 > 0 && a1 < c)

{

// a1 is not at bounds

new_b = e1 + y1 * (a1 - alph1) * k11 + y2 * (a2 - alph2) * k12 + bias;

}

else

{

if (a2 > 0 && a2 < c)

{

// a1 is at bounds but a2 is not.

new_b = e2 + y1 * (a1 - alph1) * k12 + y2 * (a2 - alph2) * k22 + bias;

}

else

{

// Both new Lagrange multipliers are at bound. SMO algorithm

// chooses the threshold to be halfway in between b1 and b2.

b1 = e1 + y1 * (a1 - alph1) * k11 + y2 * (a2 - alph2) * k12 + bias;

b2 = e2 + y1 * (a1 - alph1) * k12 + y2 * (a2 - alph2) * k22 + bias;

new_b = (b1 + b2) / 2;

}

delta_b = new_b - bias;

bias = new_b;

// Update error cache using new Lagrange multipliers

double t1 = y1 * (a1 - alph1);

double t2 = y2 * (a2 - alph2);

for (int i = 0; i < inputs.Length; i++)

{

if (0 < alpha[i] && alpha[i] < c)

{

double[] point = inputs[i];

errors[i] +=

t1 * kernel.Function(p1, point) +

t2 * kernel.Function(p2, point) -

delta_b;

}

errors[i1] = 0f;

errors[i2] = 0f;

// Update lagrange multipliers

alpha[i1] = a1;

alpha[i2] = a2;

return true;

}

/// <summary>

/// Computes the SVM output for a given point.

/// </summary>

///

private double compute(double[] point)

{

double sum = -bias;

for (int i = 0; i < inputs.Length; i++)

{

if (alpha[i] > 0)

sum += alpha[i] * outputs[i] * kernel.Function(inputs[i], point);

}

return sum;

}

private double computeComplexity()

{

// Compute initial value for C as the number of examples

// divided by the trace of the input sample kernel matrix.

double sum = 0.0;

for (int i = 0; i < inputs.Length; i++)

sum += kernel.Function(inputs[i], inputs[i]);

return inputs.Length / sum;

}

Using the code

In the following example, we will be training a Polynomial Kernel Support Vector Machine to recognize the XOR classification problem. The XOR function is classic example of a pattern classification problem that is not linearly separable.

    double[][] inputs =
    {
        new double[] { -1, -1},
        new double[] { -1,  1},
        new double[] {  1, -1},
        new double[] {  1,  1}
    };

    int[] xor = { -1, 1, 1, -1 };

double[][] inputs =

{

new double[] { -1, -1},

new double[] { -1, 1},

new double[] { 1, -1},

new double[] { 1, 1}

};

int[] xor = { -1, 1, 1, -1 };

Here, remember that the SVM is a margin classifier that classifies instances as either 1 or –1. So the training and expected output for the classification task should also be in this range. There are no such requirements for the inputs, though.

To create the Kernel Support Vector Machine with a Polynomial Kernel, do:

    // Create Kernel Support Vector Machine with a Polynomial Kernel of 2nd degree
    KernelSupportVectorMachine machine = new KernelSupportVectorMachine(
        new Polynomial(2), inputs.Length);

// Create Kernel Support Vector Machine with a Polynomial Kernel of 2nd degree

KernelSupportVectorMachine machine = new KernelSupportVectorMachine(

new Polynomial(2), inputs.Length);

After the machine has been created, create a new Learning algorithm. As we are going to do classification, we will be using the standard SequentialMinimalOptimization algorithm.

// Create the sequential minimal optimization teacher
var learn = new SequentialMinimalOptimization(machine, inputs, xor);

// Run the learning algorithm
double error = learn.Run();

// Create the sequential minimal optimization teacher

var learn = new SequentialMinimalOptimization(machine, inputs, xor);

// Run the learning algorithm

double error = learn.Run();

After the model has been trained, we can compute its outputs for the given inputs.

double[] output = machine.Compute(inputs);

1	double[] output = machine.Compute(inputs);

The machine should be able to correctly identify all of the input instances.

Sample application

The sample application is able to perform both Classification and Regression using Support Vector Machines. It can read Excel spreadsheets and determines the task to be performed depending on the number of the columns in the sheet. If the input table contains two columns (e.g. X and Y) it will be interpreted as a regression problem X –> Y. If the input table contains three columns (e.g. x1, x2 and Y) it will be interpreted as a classification problem <x1,x2> belongs to class Y, Y being either 1 or -1.

Classification

To perform classification, load a classification task data such as the Yin Yang classification problem.

Yin Yang classification problem. The goal is to create a model which best determines whether a given point belongs to class blue or green. It is a clear example of a non-linearly separable problem.

Creation of a Gaussian Kernel Support Vector Machine with σ = 1.2236, C = 1.0, ε = 0.001 and T = 0.001.

Classification using the created Support Vector Machine. Notice it achieves an accuracy of 97%, with sensitivity and specifity rates of 98% and 96%, respectively.

Regression

To perform regression, we can load the Gaussian noise sine wave example.

Noise sine wave regression problem.

Creation of a Gaussian Kernel Support Vector Machine with σ = 1.2236, C = 1.0, ε = 0.2 and T = 0.001.

After the model has been created, we can plot the model approximation for the sine wave data. The blue line shows the curve approximation for the original red training dots.

Regression using the created Kernel Support Vector Machine. Notice the coefficient of determination r² of 0.95. The closer to one, the better.

References

Wikipedia contributors, “Sequential Minimal Optimization,” Wikipedia, The Free Encyclopedia, http://en.wikipedia.org/wiki/Sequential_Minimal_Optimization (accessed April 24, 2010).
Wikipedia contributors, “Support Vector Machine,” Wikipedia, The Free Encyclopedia, http://en.wikipedia.org/wiki/Support_vector_machine (accessed April 24, 2010).
John C. Platt, Sequential Minimal Optimization: A Fast Algorithm for Training Support Vector Machines, Microsoft Research, 1998.
J. P. Lewis, A Short SVM (Support Vector Machine) Tutorial. CGIT Lab / IMSC, University of Southern California.
A. J. Smola and B. Schoelkopf. A Tutorial on Support Vector Regression. NeuroCOLT2 Technical Report Series, 1998.
S.K. Shevade et al. Improvements to SMO Algorithm for SVM Regression, 1999.
G. W. Flake, S. Lawrence. Efficient SVM Regression Training with SMO.

Partial Least Squares Analysis and Regression in C#

April 10, 2010 / cesarsouza / 17 Comments

Partial Least Squares Regression (PLS) is a technique that generalizes and combines features from principal component analysis and (multivariate) multiple regression. It has been widely adopted in the field of chemometrics and social sciences.

Download sample application and source code
Browse the sample application source code online
Browse the Partial Least Squares source code online
Download the Accord.NET Statistics Framework

The code presented here is also part of the Accord.NET Framework. The Accord.NET Framework is a framework for developing machine learning, computer vision, computer audition, statistics and math applications. Please see the starting guide for mode details. The latest version of the framework includes the latest version of this code plus many other statistics and machine learning tools.

Introduction
Overview
Source Code
Using the code
Sample application
See also
References

Introduction

Partial least squares regression (PLS-regression) is a statistical method that is related to principal components regression. The goal of this method is to find a linear regression model by projecting both the predicted variables and the observable variables to new, latent variable spaces. It was developed in the 1960s by Herman Wold to be used in econometrics. However, today it is most commonly used for regression in the field of chemometrics.

In statistics, latent variables (as opposed to observable variables), are variables
that are not directly observed but are rather inferred (through a mathematical model) from other variables that are observed (directly measured). Mathematical models that aim to explain observed variables in terms of latent variables are called latent variable models.

A PLS model will try to find the multidimensional direction in the X space that explains the maximum multidimensional variance direction in the Y space. PLS-regression is particularly suited when the matrix of predictors has more variables
than observations, and when there is multicollinearity among X values. Its interesting to note that standard linear regression would likely fail to produce meaningful interpretable models in those cases.

Overview

Multivariate Linear Regression in Latent Space

Multiple Linear Regression is a generalization of simple linear regression for multiple inputs. In turn, Multivariate Linear Regression is a generalization of Multiple Linear Regression for multiple outputs. The multivariate linear regression is a general linear regression model which can map an arbitrary dimension space into another arbitrary dimension space using only linear relationships. In the context of PLS, it is used to map the latent variable space for the inputs X into the latent variable space for the output variables Y. Those latent variable spaces are spawned by the loading matrices for X and Y, commonly denoted P and Q, respectively.

The goal of PLS algorithms are therefore to find those two matrices. There are mainly two algorithms to do this: NIPALS and SIMPLS.

Algorithm 1: NIPALS

Here is an exposition of the NIPALS algorithm for finding the loading matrices required for PLS regression. There are, however, many variations of this algorithm which normalize or do not normalize certain vectors.

Algorithm:

Let X be the mean-centered input matrix,

Let Y be the mean-centered output matrix,

Let P be the loadings matrix for X, and let p_idenote the i-th column of P;

Let Q be the loadings matrix for Y, and let q_idenote the i-th column of Q;

Let T be the score matrix for X, and t_idenote the i-th column of T;

Let U be the score matrix for Y, and u_idenote the i-th column of U;

Let W be the PLS weight matrix, and w_idenote the i-th column of W; and

Let B be a diagonal matrix of diagonal coefficients b_i

Then:

For each factor i to be calculated:

Initially choose u_i as the largest column vector in
X (having the largest sum of squares)

While (t_i has not converged to a desired precision)

w_i ∝ X’u_i(estimate X weights)

t_i ∝ Xw_i(estimate X factor scores)

q_i ∝ Y’t_i(estimate Y weights)

u_i = Yq_i(estimate Y scores)

b_i = t’u (compute prediction coefficient b)

p_i = X’t (estimate X factor loadings)

X = X – tp’ (deflate X)

In other statistical analysis such as PCA, it is often interesting to inspect how much of the variance can be explained by each of the principal component dimensions. The same can also be accomplished for PLS, both for the input (predictor) variables X and outputs (regressor) variables Y. For the input variables, the amount of variance explained by each factor can be computed as b_i². For outputs, it can be computed as the sum of the squared elements of its column in the matrix P, i.e. as Sum(p_i²).

Algorithm 2: SIMPLS

SIMPLS is an alternative algorithm for finding the PLS matrices P and Q that has been derived considering the true objective of maximizing the covariance between the latent factors and the output vectors. Both NIPALS and SIMPLS are equivalent when there is just one output variable Y to be regressed. However, their answers can differ when Y is multi-dimensional. However, because the construction of the weight vectors used by SIMPLS is based on the empirical variance–covariance matrix of the joint input and output variables, outliers present in the data can adversely impact its performance.

Algorithm:

Let X be the mean-centered input matrix,

Let Y be the mean-centered output matrix,

Let P be the loadings matrix for X, and let p_i denote the i-th column of P;

Let C be the loadings matrix for Y, and let c_i denote the i-th column of C;

Let T be the score matrix for X, and t_i denote the i-th column of T;

Let U be the score matrix for Y, and u_i denote the i-th column of U; and

Let W be the PLS weight matrix, and w_i denote the i-th column of W.

Then:

Create the covariance matrix C = X’Y

For each factor i to be calculated:

Perform SVD on the covariance matrix and store the first left singular vector in w_i and the first right singular value times the singular values in c_i.

t_i ∝ X*w_i   (estimate X factor scores)

p_i = X’*t_i           (estimate X factor loadings)

c_i = c_i/norm(t_i)     (estimate Y weights)

w_i = w_i/norm(t_i)     (estimate X weights)

u_i = Y*c_i   (estimate Y scores)

v_i = p_i              (form the basis vector v_i)

Make v orthogonal to the previous loadings V

Make u orthogonal to the previous scores T

Deflate the covariance matrix C

C = C – v_i*(v_i‘*C)

Source Code

This section contains the realization of the NIPALS and SIMPLS algorithms in C#. The models have been implemented considering an object-oriented structure that is particularly suitable to be data-bound to Windows.Forms (or WPF) controls.

Class Diagram

Class diagram for the Partial Least Squares Analysis.

Performing PLS using NIPALS

Here is the source code for computing PLS using the NIPALS algorithm:

    /// <summary>
    ///   Computes PLS parameters using NIPALS algorithm.
    /// </summary>
    /// 
    private void nipals(double[,] X, double[,] Y, int factors, double tolerance)
    {
        // References:
        //  - Hervé Abdi, http://www.utdallas.edu/~herve/abdi-wireCS-PLS2010.pdf

        // Initialize and prepare the data
        int rows = sourceX.GetLength(0);
        int xcols = sourceX.GetLength(1);
        int ycols = sourceY.GetLength(1);

        // Initialize storage variables
        double[,] E = (double[,])X.Clone();
        double[,] F = (double[,])Y.Clone();

        double[,] T = new double[rows, factors];
        double[,] U = new double[rows, factors];
        double[,] P = new double[xcols, factors];
        double[,] C = new double[ycols, factors];
        double[,] W = new double[xcols, xcols];
        double[] B = new double[xcols];

        double[] varX = new double[factors];
        double[] varY = new double[factors];

        // Initialize the algorithm
        bool stop = false;

        #region NIPALS
        for (int factor = 0; factor < factors && !stop; factor++)
        {
            // Select t as the largest column from X,
            double[] t = E.GetColumn(largest(E));

            // Select u as the largest column from Y.
            double[] u = F.GetColumn(largest(F));

            // Will store weights for X and Y
            double[] w = new double[xcols];
            double[] c = new double[ycols];

            double norm_t = Norm.Euclidean(t);

            #region Iteration
            while (norm_t > 1e-14)
            {
                // Store initial t to check convergence
                double[] t0 = (double[])t.Clone();


                // Step 1. Estimate w (X weights): w ∝ E'*u
                //   (in Abdi's paper, X is referred as E).

                // 1.1. Compute w = E'*u;
                w = new double[xcols];
                for (int j = 0; j < w.Length; j++)
                    for (int i = 0; i < u.Length; i++)
                        w[j] += E[i, j] * u[i];

                // 1.2. Normalize w (w = w/norm(w))
                w = w.Divide(Norm.Euclidean(w));


                // Step 2. Estimate t (X factor scores): t ∝ E*w
                //   (in Abdi's paper, X is referred as E).

                // 2.1. Compute t = E*w
                t = new double[rows];
                for (int i = 0; i < t.Length; i++)
                    for (int j = 0; j < w.Length; j++)
                        t[i] += E[i, j] * w[j];

                // 2.2. Normalize t: t = t/norm(t)
                t = t.Divide(norm_t = Norm.Euclidean(t));


                // Step 3. Estimate c (Y weights): c ∝ F't
                //   (in Abdi's paper, Y is referred as F).

                // 3.1. Compute c = F'*t0;
                c = new double[ycols];
                for (int j = 0; j < c.Length; j++)
                    for (int i = 0; i < t.Length; i++)
                        c[j] += F[i, j] * t[i];

                // 3.2. Normalize q: c = c/norm(q)
                c = c.Divide(Norm.Euclidean(c));


                // Step 4. Estimate u (Y scores): u = F*q
                //   (in Abdi's paper, Y is referred as F).

                // 4.1. Compute u = F*q;
                u = new double[rows];
                for (int i = 0; i < u.Length; i++)
                    for (int j = 0; j < c.Length; j++)
                        u[i] += F[i, j] * c[j];


                // Recalculate norm of the difference
                norm_t = 0.0;
                for (int i = 0; i < t.Length; i++)
                {
                    double d = (t0[i] - t[i]);
                    norm_t += d * d;
                }

                norm_t = Math.Sqrt(norm_t);
            }
            #endregion


            // Compute the value of b which is used to
            // predict Y from t as b = t'u [Abdi, 2010]
            double b = t.InnerProduct(u);

            // Compute factor loadings for X as p = E'*t [Abdi, 2010]
            double[] p = new double[xcols];
            for (int j = 0; j < p.Length; j++)
                for (int i = 0; i < rows; i++)
                    p[j] += E[i, j] * t[i];

            // Perform deflaction of X and Y
            for (int i = 0; i < t.Length; i++)
            {
                // Deflate X as X = X - t*p';
                for (int j = 0; j < p.Length; j++)
                    E[i, j] -= t[i] * p[j];

                // Deflate Y as Y = Y - b*t*q';
                for (int j = 0; j < c.Length; j++)
                    F[i, j] -= b * t[i] * c[j];
            }

            // Calculate explained variances
            varY[factor] = b * b;
            varX[factor] = p.InnerProduct(p);

            // Save iteration results
            T.SetColumn(factor, t);
            P.SetColumn(factor, p);
            U.SetColumn(factor, u);
            C.SetColumn(factor, c);
            W.SetColumn(factor, w);
            B[factor] = b;

            // Check for residuals as stop criteria
            double[] norm_x = Norm.Euclidean(E);
            double[] norm_y = Norm.Euclidean(F);

            stop = true;
            for (int i = 0; i < norm_x.Length && stop == true; i++)
            {
                // If any of the residuals is higher than the tolerance
                if (norm_x[i] > tolerance || norm_y[i] > tolerance)
                    stop = false;
            }
        }

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

/// <summary>

/// Computes PLS parameters using NIPALS algorithm.

/// </summary>

///

private void nipals(double[,] X, double[,] Y, int factors, double tolerance)

{

// References:

// - Hervé Abdi, http://www.utdallas.edu/~herve/abdi-wireCS-PLS2010.pdf

// Initialize and prepare the data

int rows = sourceX.GetLength(0);

int xcols = sourceX.GetLength(1);

int ycols = sourceY.GetLength(1);

// Initialize storage variables

double[,] E = (double[,])X.Clone();

double[,] F = (double[,])Y.Clone();

double[,] T = new double[rows, factors];

double[,] U = new double[rows, factors];

double[,] P = new double[xcols, factors];

double[,] C = new double[ycols, factors];

double[,] W = new double[xcols, xcols];

double[] B = new double[xcols];

double[] varX = new double[factors];

double[] varY = new double[factors];

// Initialize the algorithm

bool stop = false;

#region NIPALS

for (int factor = 0; factor < factors && !stop; factor++)

{

// Select t as the largest column from X,

double[] t = E.GetColumn(largest(E));

// Select u as the largest column from Y.

double[] u = F.GetColumn(largest(F));

// Will store weights for X and Y

double[] w = new double[xcols];

double[] c = new double[ycols];

double norm_t = Norm.Euclidean(t);

#region Iteration

while (norm_t > 1e-14)

{

// Store initial t to check convergence

double[] t0 = (double[])t.Clone();

// Step 1. Estimate w (X weights): w ∝ E'*u

// (in Abdi's paper, X is referred as E).

// 1.1. Compute w = E'*u;

w = new double[xcols];

for (int j = 0; j < w.Length; j++)

for (int i = 0; i < u.Length; i++)

w[j] += E[i, j] * u[i];

// 1.2. Normalize w (w = w/norm(w))

w = w.Divide(Norm.Euclidean(w));

// Step 2. Estimate t (X factor scores): t ∝ E*w

// (in Abdi's paper, X is referred as E).

// 2.1. Compute t = E*w

t = new double[rows];

for (int i = 0; i < t.Length; i++)

for (int j = 0; j < w.Length; j++)

t[i] += E[i, j] * w[j];

// 2.2. Normalize t: t = t/norm(t)

t = t.Divide(norm_t = Norm.Euclidean(t));

// Step 3. Estimate c (Y weights): c ∝ F't

// (in Abdi's paper, Y is referred as F).

// 3.1. Compute c = F'*t0;

c = new double[ycols];

for (int j = 0; j < c.Length; j++)

for (int i = 0; i < t.Length; i++)

c[j] += F[i, j] * t[i];

// 3.2. Normalize q: c = c/norm(q)

c = c.Divide(Norm.Euclidean(c));

// Step 4. Estimate u (Y scores): u = F*q

// (in Abdi's paper, Y is referred as F).

// 4.1. Compute u = F*q;

u = new double[rows];

for (int i = 0; i < u.Length; i++)

for (int j = 0; j < c.Length; j++)

u[i] += F[i, j] * c[j];

// Recalculate norm of the difference

norm_t = 0.0;

for (int i = 0; i < t.Length; i++)

{

double d = (t0[i] - t[i]);

norm_t += d * d;

}

norm_t = Math.Sqrt(norm_t);

}

#endregion

// Compute the value of b which is used to

// predict Y from t as b = t'u [Abdi, 2010]

double b = t.InnerProduct(u);

// Compute factor loadings for X as p = E'*t [Abdi, 2010]

double[] p = new double[xcols];

for (int j = 0; j < p.Length; j++)

for (int i = 0; i < rows; i++)

p[j] += E[i, j] * t[i];

// Perform deflaction of X and Y

for (int i = 0; i < t.Length; i++)

{

// Deflate X as X = X - t*p';

for (int j = 0; j < p.Length; j++)

E[i, j] -= t[i] * p[j];

// Deflate Y as Y = Y - b*t*q';

for (int j = 0; j < c.Length; j++)

F[i, j] -= b * t[i] * c[j];

}

// Calculate explained variances

varY[factor] = b * b;

varX[factor] = p.InnerProduct(p);

// Save iteration results

T.SetColumn(factor, t);

P.SetColumn(factor, p);

U.SetColumn(factor, u);

C.SetColumn(factor, c);

W.SetColumn(factor, w);

B[factor] = b;

// Check for residuals as stop criteria

double[] norm_x = Norm.Euclidean(E);

double[] norm_y = Norm.Euclidean(F);

stop = true;

for (int i = 0; i < norm_x.Length && stop == true; i++)

{

// If any of the residuals is higher than the tolerance

if (norm_x[i] > tolerance || norm_y[i] > tolerance)

stop = false;

}

Performing PLS using SIMPLS

And here is the source code for computing PLS using the SIMPLS algorithm:

    /// <summary>
    ///   Computes PLS parameters using SIMPLS algorithm.
    /// </summary>
    /// 
    private void simpls(double[,] X, double[,] Y, int factors)
    {
        // References:
        //  - Martin Anderson, "A comparison of nine PLS1 algorithms". Journal of Chemometrics,
        //    2009. Available on: http://onlinelibrary.wiley.com/doi/10.1002/cem.1248/pdf
        //  - Hervé Abdi, http://www.utdallas.edu/~herve/abdi-wireCS-PLS2010.pdf
        //  - Statsoft, http://www.statsoft.com/textbook/partial-least-squares/#SIMPLS
        //  - Sijmen de Jong, "SIMPLS: an alternative approach to partial least squares regression"
        //  - N.M. Faber and J. Ferré, “On the numerical stability of two widely used PLS algorithms,”
        //    J. Chemometrics, 22, pps 101-105, 2008.

        // Initialize and prepare the data
        int rows = sourceX.GetLength(0);
        int xcols = sourceX.GetLength(1);
        int ycols = sourceY.GetLength(1);

        // Initialize storage variables
        double[,] P = new double[xcols, factors]; // loading matrix P, the loadings for X such that X = TP + F
        double[,] C = new double[ycols, factors]; // loading matrix C, the loadings for Y such that Y = TC + E
        double[,] T = new double[rows, factors];  // factor score matrix T
        double[,] U = new double[rows, factors];  // factor score matrix U
        double[,] W = new double[xcols, factors]; // weight matrix W

        double[] varX = new double[factors];
        double[] varY = new double[factors];

        // Orthogonal loadings
        double[,] V = new double[xcols, factors];

        // Create covariance matrix C = X'Y
        double[,] covariance = X.TransposeAndMultiply(Y);

        #region SIMPLS
        for (int factor = 0; factor < factors; factor++)
        {
            // Step 1. Obtain the dominant eigenvector w of C'C. However, we
            //   can avoid computing the matrix multiplication by using the
            //   singular value decomposition instead, which is also more
            //   stable. the ﬁrst weight vector w is the left singular vector
            //   of C=X'Y [Abdi, 2007].

            var svd = new SingularValueDecomposition(covariance,
                computeLeftSingularVectors: true,
                computeRightSingularVectors: false,
                autoTranspose: true);

            double[] w = svd.LeftSingularVectors.GetColumn(0);
            double[] c = covariance.TransposeAndMultiply(w);


            // Step 2. Estimate X factor scores: t ∝ X*w
            //   Similarly to NIPALS, the T factor of SIMPLS
            //   is computed as T=X*W [Statsoft] [Abdi, 2010].

            // 2.1. Estimate t (X factor scores): t = X*w [Abdi, 2010]
            double[] t = new double[rows];
            for (int i = 0; i < t.Length; i++)
                for (int j = 0; j < w.Length; j++)
                    t[i] += X[i, j] * w[j];

            // 2.2. Normalize t (X factor scores): t = t/norm(t)
            double norm_t = Norm.Euclidean(t);
            t = t.Divide(norm_t);


            // Step 3. Estimate p (X factor loadings): p = X'*t
            double[] p = new double[xcols];
            for (int i = 0; i < p.Length; i++)
                for (int j = 0; j < t.Length; j++)
                    p[i] += X[j, i] * t[j];


            // Step 4. Estimate X and Y weights. Actually, the weights have
            //   been computed in the first step during SVD. However, since
            //   the X factor scores have been normalized, we also have to
            //   normalize weights accordingly: w = w/norm(t), c = c/norm(t)
            w = w.Divide(norm_t);
            c = c.Divide(norm_t);


            // Step 5. Estimate u (Y factor scores): u = Y*c [Abdi, 2010]
            double[] u = new double[rows];
            for (int i = 0; i < u.Length; i++)
                for (int j = 0; j < c.Length; j++)
                    u[i] += Y[i, j] * c[j];


            // Step 6. Initialize the orthogonal loadings
            double[] v = (double[])p.Clone();


            // Step 7. Make v orthogonal to the previous loadings
            // http://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process

            if (factor > 0)
            {
                // 7.1. MGS for v [Martin Anderson, 2009]
                for (int j = 0; j < factor; j++)
                {
                    double proj = 0.0;
                    for (int k = 0; k < v.Length; k++)
                        proj += v[k] * V[k, j];

                    for (int k = 0; k < v.Length; k++)
                        v[k] -= proj * V[k, j];
                }

                // 7.1. MGS for u [Martin Anderson, 2009]
                for (int j = 0; j < factor; j++)
                {
                    double proj = 0.0;
                    for (int k = 0; k < u.Length; k++)
                        proj += u[k] * T[k, j];

                    for (int k = 0; k < u.Length; k++)
                        u[k] -= proj * T[k, j];
                }
            }

            // 7.2. Normalize orthogonal loadings
            v = v.Divide(Norm.Euclidean(v));


            // Step 8. Deflate covariance matrix as s = s - v * (v' * s)
            //   as shown in simpls1 in [Martin Anderson, 2009] appendix.
            double[,] cov = (double[,])covariance.Clone();
            for (int i = 0; i < v.Length; i++)
            {
                for (int j = 0; j < v.Length; j++)
                {
                    double d = v[i] * v[j];

                    for (int k = 0; k < ycols; k++)
                        cov[i, k] -= d * covariance[j, k];
                }
            }
            covariance = cov;

            // Save iteration
            W.SetColumn(factor, w);
            U.SetColumn(factor, u);
            C.SetColumn(factor, c);
            T.SetColumn(factor, t);
            P.SetColumn(factor, p);
            V.SetColumn(factor, v);

            // Compute explained variance
            varX[factor] = p.InnerProduct(p);
            varY[factor] = c.InnerProduct(c);
        }
        #endregion

        // Set class variables
        this.scoresX = T;      // factor score matrix T
        this.scoresY = U;      // factor score matrix U
        this.loadingsX = P;    // loading matrix P, the loadings for X such that X = TP + F
        this.loadingsY = C;    // loading matrix Q, the loadings for Y such that Y = TQ + E
        this.weights = W;      // the columns of R are weight vectors
        this.coeffbase = W;

        // Calculate variance explained proportions
        this.componentProportionX = new double[factors];
        this.componentProportionY = new double[factors];

        double sumX = 0.0, sumY = 0.0;
        for (int i = 0; i < rows; i++)
        {
            // Sum of squares for matrix X
            for (int j = 0; j < xcols; j++)
                sumX += X[i, j] * X[i, j];

            // Sum of squares for matrix Y
            for (int j = 0; j < ycols; j++)
                sumY += Y[i, j] * Y[i, j];
        }

        // Calculate variance proportions
        for (int i = 0; i < factors; i++)
        {
            componentProportionY[i] = varY[i] / sumY;
            componentProportionX[i] = varX[i] / sumX;
        }
    }

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

/// <summary>

/// Computes PLS parameters using SIMPLS algorithm.

/// </summary>

///

private void simpls(double[,] X, double[,] Y, int factors)

{

// References:

// - Martin Anderson, "A comparison of nine PLS1 algorithms". Journal of Chemometrics,

// 2009. Available on: http://onlinelibrary.wiley.com/doi/10.1002/cem.1248/pdf

// - Hervé Abdi, http://www.utdallas.edu/~herve/abdi-wireCS-PLS2010.pdf

// - Statsoft, http://www.statsoft.com/textbook/partial-least-squares/#SIMPLS

// - Sijmen de Jong, "SIMPLS: an alternative approach to partial least squares regression"

// - N.M. Faber and J. Ferré, “On the numerical stability of two widely used PLS algorithms,”

// J. Chemometrics, 22, pps 101-105, 2008.

// Initialize and prepare the data

int rows = sourceX.GetLength(0);

int xcols = sourceX.GetLength(1);

int ycols = sourceY.GetLength(1);

// Initialize storage variables

double[,] P = new double[xcols, factors]; // loading matrix P, the loadings for X such that X = TP + F

double[,] C = new double[ycols, factors]; // loading matrix C, the loadings for Y such that Y = TC + E

double[,] T = new double[rows, factors]; // factor score matrix T

double[,] U = new double[rows, factors]; // factor score matrix U

double[,] W = new double[xcols, factors]; // weight matrix W

double[] varX = new double[factors];

double[] varY = new double[factors];

// Orthogonal loadings

double[,] V = new double[xcols, factors];

// Create covariance matrix C = X'Y

double[,] covariance = X.TransposeAndMultiply(Y);

#region SIMPLS

for (int factor = 0; factor < factors; factor++)

{

// Step 1. Obtain the dominant eigenvector w of C'C. However, we

// can avoid computing the matrix multiplication by using the

// singular value decomposition instead, which is also more

// stable. the ﬁrst weight vector w is the left singular vector

// of C=X'Y [Abdi, 2007].

var svd = new SingularValueDecomposition(covariance,

computeLeftSingularVectors: true,

computeRightSingularVectors: false,

autoTranspose: true);

double[] w = svd.LeftSingularVectors.GetColumn(0);

double[] c = covariance.TransposeAndMultiply(w);

// Step 2. Estimate X factor scores: t ∝ X*w

// Similarly to NIPALS, the T factor of SIMPLS

// is computed as T=X*W [Statsoft] [Abdi, 2010].

// 2.1. Estimate t (X factor scores): t = X*w [Abdi, 2010]

double[] t = new double[rows];

for (int i = 0; i < t.Length; i++)

for (int j = 0; j < w.Length; j++)

t[i] += X[i, j] * w[j];

// 2.2. Normalize t (X factor scores): t = t/norm(t)

double norm_t = Norm.Euclidean(t);

t = t.Divide(norm_t);

// Step 3. Estimate p (X factor loadings): p = X'*t

double[] p = new double[xcols];

for (int i = 0; i < p.Length; i++)

for (int j = 0; j < t.Length; j++)

p[i] += X[j, i] * t[j];

// Step 4. Estimate X and Y weights. Actually, the weights have

// been computed in the first step during SVD. However, since

// the X factor scores have been normalized, we also have to

// normalize weights accordingly: w = w/norm(t), c = c/norm(t)

w = w.Divide(norm_t);

c = c.Divide(norm_t);

// Step 5. Estimate u (Y factor scores): u = Y*c [Abdi, 2010]

double[] u = new double[rows];

for (int i = 0; i < u.Length; i++)

for (int j = 0; j < c.Length; j++)

u[i] += Y[i, j] * c[j];

// Step 6. Initialize the orthogonal loadings

double[] v = (double[])p.Clone();

// Step 7. Make v orthogonal to the previous loadings

// http://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process

if (factor > 0)

{

// 7.1. MGS for v [Martin Anderson, 2009]

for (int j = 0; j < factor; j++)

{

double proj = 0.0;

for (int k = 0; k < v.Length; k++)

proj += v[k] * V[k, j];

for (int k = 0; k < v.Length; k++)

v[k] -= proj * V[k, j];

}

// 7.1. MGS for u [Martin Anderson, 2009]

for (int j = 0; j < factor; j++)

{

double proj = 0.0;

for (int k = 0; k < u.Length; k++)

proj += u[k] * T[k, j];

for (int k = 0; k < u.Length; k++)

u[k] -= proj * T[k, j];

}

// 7.2. Normalize orthogonal loadings

v = v.Divide(Norm.Euclidean(v));

// Step 8. Deflate covariance matrix as s = s - v * (v' * s)

// as shown in simpls1 in [Martin Anderson, 2009] appendix.

double[,] cov = (double[,])covariance.Clone();

for (int i = 0; i < v.Length; i++)

{

for (int j = 0; j < v.Length; j++)

{

double d = v[i] * v[j];

for (int k = 0; k < ycols; k++)

cov[i, k] -= d * covariance[j, k];

}

covariance = cov;

// Save iteration

W.SetColumn(factor, w);

U.SetColumn(factor, u);

C.SetColumn(factor, c);

T.SetColumn(factor, t);

P.SetColumn(factor, p);

V.SetColumn(factor, v);

// Compute explained variance

varX[factor] = p.InnerProduct(p);

varY[factor] = c.InnerProduct(c);

}

#endregion

// Set class variables

this.scoresX = T; // factor score matrix T

this.scoresY = U; // factor score matrix U

this.loadingsX = P; // loading matrix P, the loadings for X such that X = TP + F

this.loadingsY = C; // loading matrix Q, the loadings for Y such that Y = TQ + E

this.weights = W; // the columns of R are weight vectors

this.coeffbase = W;

// Calculate variance explained proportions

this.componentProportionX = new double[factors];

this.componentProportionY = new double[factors];

double sumX = 0.0, sumY = 0.0;

for (int i = 0; i < rows; i++)

{

// Sum of squares for matrix X

for (int j = 0; j < xcols; j++)

sumX += X[i, j] * X[i, j];

// Sum of squares for matrix Y

for (int j = 0; j < ycols; j++)

sumY += Y[i, j] * Y[i, j];

}

// Calculate variance proportions

for (int i = 0; i < factors; i++)

{

componentProportionY[i] = varY[i] / sumY;

componentProportionX[i] = varX[i] / sumX;

}

Multivariate Linear Regression

Multivariate Linear Regression is computed in a similar manner to a Multiple Linear Regression. The only difference is that, instead of having a weight vector and a intercept, we have a weight matrix and a intercept vector.

    /// <summary> 
    /// Computes the Multiple Linear Regression for an input vector. 
    /// </summary> 
    /// <param name="input">The input vector.</param> 
    /// <returns>The calculated output.</returns> 
    public double[] Compute(double[] input) 
    { 
        int N = input.Length; 
        int M = coefficients.GetLength(1); 

        double[] result = new double[M]; 
        for (int i = 0; i < M; i++) 
        { 
            result[i] = intercepts[i]; 

            for (int j = 0; j < N; j++) 
                result[i] += input[j] * coefficients[j, i]; 
        } 

        return result; 
    }

/// <summary>

/// Computes the Multiple Linear Regression for an input vector.

/// </summary>

/// <param name="input">The input vector.</param>

/// <returns>The calculated output.</returns>

public double[] Compute(double[] input)

{

int N = input.Length;

int M = coefficients.GetLength(1);

double[] result = new double[M];

for (int i = 0; i < M; i++)

{

result[i] = intercepts[i];

for (int j = 0; j < N; j++)

result[i] += input[j] * coefficients[j, i];

}

return result;

}

The weight matrix and the intercept vector are computed in the PartialLeastSquaresAnalysis class by the CreateRegression method. In case the analyzed data already was mean centered before being fed to the analysis, the constructed intercept vector will consist only of zeros.

    /// <summary> 
    /// Creates a Multivariate Linear Regression model using 
    /// coefficients obtained by the Partial Least Squares. 
    /// </summary> 
    public MultivariateLinearRegression CreateRegression(int factors) 
    { 
        int rows = sourceX.GetLength(0); 
        int xcols = sourceX.GetLength(1); 
        int ycols = sourceY.GetLength(1); 

        // Compute regression coefficients B of Y on X as B = RQ' 
        double[,] B = new double[xcols, ycols]; 
        for (int i = 0; i < xcols; i++) 
            for (int j = 0; j < ycols; j++) 
                for (int k = 0; k < factors; k++) 
                    B[i, j] += coeffbase[i, k] * loadingsY[j, k]; 

        // Divide by standard deviation if X has been normalized 
        if (analysisMethod == AnalysisMethod.Correlation) 
            for (int i = 0; i < xcols; i++) 
                for (int j = 0; j < ycols; j++) 
                    B[i, j] = B[i, j] / stdDevX[i]; 

        // Compute regression intercepts A as A = meanY - meanX'*B 
        double[] A = new double[ycols]; 
        for (int i = 0; i < ycols; i++) 
         { 
            double sum = 0.0; 
            for (int j = 0; j < xcols; j++) 
                sum += meanX[j] * B[j, i]; 
            A[i] = meanY[i] - sum; 
        } 

        return new MultivariateLinearRegression(B, A, true); 
    }

/// <summary>

/// Creates a Multivariate Linear Regression model using

/// coefficients obtained by the Partial Least Squares.

/// </summary>

public MultivariateLinearRegression CreateRegression(int factors)

{

int rows = sourceX.GetLength(0);

int xcols = sourceX.GetLength(1);

int ycols = sourceY.GetLength(1);

// Compute regression coefficients B of Y on X as B = RQ'

double[,] B = new double[xcols, ycols];

for (int i = 0; i < xcols; i++)

for (int j = 0; j < ycols; j++)

for (int k = 0; k < factors; k++)

B[i, j] += coeffbase[i, k] * loadingsY[j, k];

// Divide by standard deviation if X has been normalized

if (analysisMethod == AnalysisMethod.Correlation)

for (int i = 0; i < xcols; i++)

for (int j = 0; j < ycols; j++)

B[i, j] = B[i, j] / stdDevX[i];

// Compute regression intercepts A as A = meanY - meanX'*B

double[] A = new double[ycols];

for (int i = 0; i < ycols; i++)

{

double sum = 0.0;

for (int j = 0; j < xcols; j++)

sum += meanX[j] * B[j, i];

A[i] = meanY[i] - sum;

}

return new MultivariateLinearRegression(B, A, true);

}

Using the code

As an example, lets consider the example data from Hervé Abdi, where the goal is to predict the subjective evaluation of a set of 5 wines. The dependent variables that we want to predict for each wine are its likeability, and how well it goes with meat, or dessert (as rated by a panel of experts). The predictors are the price, the sugar, alcohol, and acidity content of each wine.

    double[,] dependent = 
    { 
        // Wine | Hedonic | Goes with meat | Goes with dessert 
        {           14,          7,                 8 }, 
        {           10,          7,                 6 }, 
        {            8,          5,                 5 }, 
        {            2,          4,                 7 }, 
        {            6,          2,                 4 }, 
    }; 

    double[,] predictors = 
    { 
        // Wine | Price | Sugar | Alcohol | Acidity 
        {   7,     7,      13,        7 }, 
        {   4,     3,      14,        7 }, 
        {  10,     5,      12,        5 }, 
        {  16,     7,      11,        3 }, 
        {  13,     3,      10,        3 }, 
    };

double[,] dependent =

{

// Wine | Hedonic | Goes with meat | Goes with dessert

{ 14, 7, 8 },

{ 10, 7, 6 },

{ 8, 5, 5 },

{ 2, 4, 7 },

{ 6, 2, 4 },

};

double[,] predictors =

{

// Wine | Price | Sugar | Alcohol | Acidity

{ 7, 7, 13, 7 },

{ 4, 3, 14, 7 },

{ 10, 5, 12, 5 },

{ 16, 7, 11, 3 },

{ 13, 3, 10, 3 },

};

Next, we proceed to create the Partial Least Squares Analysis using the Covariance
method (data will only be mean centered but not normalized) and using the SIMPLS
algorithm.

// Create the analysis using the covariance method and the SIMPLS algorithm 
PartialLeastSquaresAnalysis pls = new PartialLeastSquaresAnalysis(predictors, dependent, 
    AnalysisMethod.Covariance, PartialLeastSquaresAlgorithm.SIMPLS); 

// Compute the analysis 
pls.Compute();

// Create the analysis using the covariance method and the SIMPLS algorithm

PartialLeastSquaresAnalysis pls = new PartialLeastSquaresAnalysis(predictors, dependent,

AnalysisMethod.Covariance, PartialLeastSquaresAlgorithm.SIMPLS);

// Compute the analysis

pls.Compute();

After the analysis has been computed, we can proceed and create the regression model.

// Create the Multivariate Linear Regression model
MultivariateLinearRegression regression = pls.CreateRegression();

// Compute the regression outputs for the predictor variables
double[][] aY = regression.Compute(predictors.ToArray());

// Create the Multivariate Linear Regression model

MultivariateLinearRegression regression = pls.CreateRegression();

// Compute the regression outputs for the predictor variables

double[][] aY = regression.Compute(predictors.ToArray());

Now after the regression has been computed, we can find how well it has performed. The coefficient of determination r² for the variables Hedonic, Goes with Meat and Goes with Dessert can be computed by the CoefficientOfDetermination method of the MultivariateRegressionClass and will be, respectively, 0.9999, 0.9999 and 0.8750 – the closer to one, the better.

Sample application

The accompanying sample application performs Partial Least Squares Analysis and Regression in Excel worksheets. The predictors and dependent variables can be selected once the data has been loaded in the application.

Left: Wine example from Hervé Abdi. Right: Variance explained by PLS using the SIMPLS algorithm

Left: Partial Least Squares Analysis results and regression coefficients for the
full regression model. Right: projection of the dependent and predictors variables
using the three first factors.

Results from the Multivariate Linear Regression performed in Latent Space using
three factors from PLS.

Left: Example data from Geladi and Kowalski. Right: Analysis results for the data using NIPALS. We can see that just two factors are enough to explain the whole variance of the set.

Left: Projections to the latent spaces. Right: Loadings and coefficients for the
two factor Multivariate Linear Regression model.

Results from the Partial Least Squares Regression showing a perfect fit using only the first two factors.

References

Abdi, H. Partial least square regression, projection on latent structure regression, PLS-Regression. Wiley Interdisciplinary Reviews: Computational Statistics, 2, 97-106, 2010.
Abdi, H. Partial least square regression (PLS regression). In N.J. Salkind (Ed.): Encyclopedia of Measurement and Statistics. Thousand Oaks (CA): Sage. pp. 740-744, 2007.
Mevik, B-H. and Wehrens, R. The pls Package: Principal Component and Partial Least Squares Regression in R. Journal of Statistical Software, Volume 18, Issue 2, 2007.
De Jong, S. SIMPLS: an alternative approach to partial least squares regression. Chemometrics and Intelligent Laboratory Systems, 18: 251–263, 1993.
Geladi, P. and Kowalski, B.R.
An example of 2-block predictive partial least squares regression with simulated
data. Anal. Chim. Acta 185, 19-32 (1986).
Rosipal, R. and Trejo, L.J.
Kernel Partial Least Squares for Nonlinear Regression and Discrimination.
The Journal of Machine Learning Research, Volume 2, pp. 123, 2002.
Garson, D. Partial Least Squares Regression (PLS). Website.
Wikipedia contributors. “Partial least squares.” Wikipedia, The Free Encyclopedia. Website.

Hidden Markov Model -Based Sequence Classifiers in C#

March 24, 2010 / cesarsouza / 58 Comments

Distinct Hidden Markov Models can be trained individually on data obtained from each class of a classification problem. If we create a set of those models and specialize each model to recognize each of the separate classes, we can then use the HMMs ability to calculate the likelihood that a given sequence belongs to the model to determine the most likely class of an unknown sequence.

Download sample and source code
Download the Accord.NET Framework

This post is a continuation for the previous
article, Hidden Markov Models in C#.

The code presented here is also part of the Accord.NET Framework. The Accord.NET Framework is a framework for developing machine learning, computer vision, computer audition, statistics and math applications. It is based on the already excellent AForge.NET Framework. Please see the starting guide for mode details. The latest version of the framework includes the latest version of this code plus many other statistics and machine learning tools.

Introduction

Hidden Markov Models

An introduction about Hidden Markov Models (HMM) was presented in a previous article, entitled Hidden Markov Models in C#. HMMs are models are stochastic methods to model temporal and sequence data.

They allow, among other things, (1) to infer the most likely sequence of states that produced a given output sequence, to (2) infer which will be the most likely next state (and thus predicting the next output) and (3) calculate the probability that a given sequence of outputs originated from the system (allowing the use of hidden Markov models for sequence classification). In the context of this article, we will be more interested in results from ability (3).

Hidden Markov Models can be seem as finite state machines where for each sequence unit observation there is a state transition and, for each state, there is a output symbol emission. The picture on the left summarizes the overall definition of a HMM given in the previous article.

Hidden Markov Model -Based Sequence Classifier

Hidden Markov Models can be trained individually on data obtained from each class of a classification problem. If we create a set of models and specialize each model to recognize each of the separated classes, then we will be able to explore the HMMs ability to calculate the likelihood that a given sequence belongs to itself to perform classification of discrete sequences.

After all models have been trained, the probability of the unknown-class sequence can be computed for each model. As each model specialized in a given class, the one which outputs the highest probability can be used to determine the most likely class for the new sequence, as shown in the picture below.

Schematic representation of the sequence classification procedure.

Source Code

The implementation is very straightforward, as shown in the picture below. The Tag property of the HiddenMarkovModel class can be used to store additional information about which class the model will be representing.

Class diagram for the HMM-based sequence classifier.

Computing the most likely class for a given sequence

Computing the most likely class for a given sequence is as simple as iterating through all models in the set and selecting the one which has produced the highest likelihood value.

<span>/// &lt;summary&gt;</span><br /><span>///   Computes the most likely class for a given sequence.</span><br /><span>/// &lt;/summary&gt;</span><br /><span>public</span> <span>int</span> Compute(<span>int</span>[] sequence, <span>out</span> <span>double</span> likelihood)<br />{<br />    <span>int</span> label = 0;<br />    likelihood = 0.0;<br /><br />    <span>// For every model in the set,</span><br />    <span>for</span> (<span>int</span> i = 0; i &lt; models.Length; i++)<br />    {<br />        <span>// Evaluate the probability for the given sequence</span><br />        <span>double</span> p = models[i].Evaluate(sequence);<br /><br />        <span>// And select the one which produces the highest probability</span><br />        <span>if</span> (p &gt; likelihood)<br />        {<br />            label = i;<br />            likelihood = p;<br />        }<br />    }<br /><br />    <span>// Returns the index of the most likely model.</span><br />    <span>return</span> label;<br />}

/// <summary> /// Computes the most likely class for a given sequence. /// </summary> public int Compute(int[] sequence, out double likelihood) { int label = 0; likelihood = 0.0; // For every model in the set, for (int i = 0; i < models.Length; i++) { // Evaluate the probability for the given sequence double p = models[i].Evaluate(sequence); // And select the one which produces the highest probability if (p > likelihood) { label = i; likelihood = p; } } // Returns the index of the most likely model. return label; }

Training each model to recognize each of the output classes

To train each model to recognize each of the output classes we have to split the training data into subsets whose outputs corresponds to each of the classes from the classification problem.

<span>/// &lt;summary&gt;</span><br /><span>///   Trains each model to recognize each of the output labels.</span><br /><span>/// &lt;/summary&gt;</span><br /><span>/// &lt;returns&gt;The sum likelihood for all models after training.&lt;/returns&gt;</span><br /><span>public</span> <span>double</span> Learn(<span>int</span>[][] inputs, <span>int</span>[] outputs, <span>int</span> iterations, <span>double</span> limit)<br />{<br />    <span>double</span> sum = 0;<br /><br />    <span>// For each model,</span><br />    <span>for</span> (<span>int</span> i = 0; i &lt; models.Length; i++)<br />    {<br />        <span>// Select the input/output set corresponding</span><br />        <span>//  to the model's specialization class</span><br />        <span>int</span>[] inx = outputs.Find(y =&gt; y == i);<br />        <span>int</span>[][] observations = inputs.Submatrix(inx);<br /><br />        <span>// Train the current model in the input/output subset</span><br />        sum += models[i].Learn(observations, iterations, limit);<br />    }<br /><br />    <span>// Returns the sum likelihood for all models.</span><br />    <span>return</span> sum;<br />}

/// <summary> /// Trains each model to recognize each of the output labels. /// </summary> /// <returns>The sum likelihood for all models after training.</returns> public double Learn(int[][] inputs, int[] outputs, int iterations, double limit) { double sum = 0; // For each model, for (int i = 0; i < models.Length; i++) { // Select the input/output set corresponding // to the model's specialization class int[] inx = outputs.Find(y => y == i); int[][] observations = inputs.Submatrix(inx); // Train the current model in the input/output subset sum += models[i].Learn(observations, iterations, limit); } // Returns the sum likelihood for all models. return sum; }

Using the Code

Code usage is very simple. Once one has the input data available in the form of a sequence array and output data in form of a integer array, create a new HiddenMarkovClassifier object with the appropriate parameters.

Then, just call Learn() passing the input and output data and the convergence limit for the learning algorithm. After that, call Compute() passing a new integer sequence to be classified by the system.

<span>// Declare some testing data</span><br /><span>int</span>[][] inputs = <span>new</span> <span>int</span>[][]<br />{<br />    <span>new</span> <span>int</span>[] { 0,1,1,0 },   <span>// Class 0</span><br />    <span>new</span> <span>int</span>[] { 0,0,1,0 },   <span>// Class 0</span><br />    <span>new</span> <span>int</span>[] { 0,1,1,1,0 }, <span>// Class 0</span><br />    <span>new</span> <span>int</span>[] { 0,1,0 },     <span>// Class 0</span><br /><br />    <span>new</span> <span>int</span>[] { 1,0,0,1 },   <span>// Class 1</span><br />    <span>new</span> <span>int</span>[] { 1,1,0,1 },   <span>// Class 1</span><br />    <span>new</span> <span>int</span>[] { 1,0,0,0,1 }, <span>// Class 1</span><br />    <span>new</span> <span>int</span>[] { 1,0,1 },     <span>// Class 1</span><br />};<br /><br /><span>int</span>[] outputs = <span>new</span> <span>int</span>[]<br />{<br />    0,0,0,0, <span>// First four sequences are of class 0</span><br />    1,1,1,1, <span>// Last four sequences are of class 1</span><br />};<br /><br /><br /><span>// We are trying to predict two different classes</span><br /><span>int</span> classes = 2;<br /><br /><span>// Each sequence may have up to two symbols (0 or 1)</span><br /><span>int</span> symbols = 2;<br /><br /><span>// Nested models will have two states each</span><br /><span>int</span>[] states = <span>new</span> <span>int</span>[] { 2, 2 };<br /><br /><span>// Creates a new Hidden Markov Model Classifier with the given parameters</span><br />HiddenMarkovClassifier hmc = <span>new</span> HiddenMarkovClassifier(classes, symbols, states);<br /><br /><br /><span>// Will train until convergence of the average likelihood</span><br /><span>double</span> likelihood = hmc.Learn(inputs, outputs, 0.001);<br /><br /><br /><span>// Will assert the models have learned the sequences correctly.</span><br /><span>for</span> (<span>int</span> i = 0; i &lt; inputs.Length; i++)<br />{<br />    <span>int</span> expected = outputs[i];<br />    <span>int</span> actual = hmc.Compute(inputs[i], <span>out</span> likelihood);<br />    Assert.AreEqual(expected, actual);<br />}

// Declare some testing data int[][] inputs = new int[][] { new int[] { 0,1,1,0 }, // Class 0 new int[] { 0,0,1,0 }, // Class 0 new int[] { 0,1,1,1,0 }, // Class 0 new int[] { 0,1,0 }, // Class 0 new int[] { 1,0,0,1 }, // Class 1 new int[] { 1,1,0,1 }, // Class 1 new int[] { 1,0,0,0,1 }, // Class 1 new int[] { 1,0,1 }, // Class 1 }; int[] outputs = new int[] { 0,0,0,0, // First four sequences are of class 0 1,1,1,1, // Last four sequences are of class 1 }; // We are trying to predict two different classes int classes = 2; // Each sequence may have up to two symbols (0 or 1) int symbols = 2; // Nested models will have two states each int[] states = new int[] { 2, 2 }; // Creates a new Hidden Markov Model Classifier with the given parameters HiddenMarkovClassifier hmc = new HiddenMarkovClassifier(classes, symbols, states); // Will train until convergence of the average likelihood double likelihood = hmc.Learn(inputs, outputs, 0.001); // Will assert the models have learned the sequences correctly. for (int i = 0; i < inputs.Length; i++) { int expected = outputs[i]; int actual = hmc.Compute(inputs[i], out likelihood); Assert.AreEqual(expected, actual); }

Sample Application

The accompanying sample application can read Excel spreadsheets containing integer sequences and labels for those sequences. An example spreadsheet contained inside the Resources folder in the compressed archive file demonstrates the expected format for this data.

Left: Input sequences data. Right: Hidden Markov Models contained in the Sequence Classifier with initial configurations.

Left: Trained models. Right: Performance test for the Hidden Markov Model Sequence Classifier.

References

Kadous, Mohammed W., “A general architecture for supervised classification of multivariate time series”, Technical Report UNSW-CSE-TR-9806. School of Computer Science & Engineering, University of New South Wales, 1998.

Hidden Markov Models in C#

March 23, 2010 / cesarsouza / 20 Comments

Hidden Markov Models (HMM) are stochastic methods to model temporal and sequence data. They are especially known for their application in temporal pattern recognition such as speech, handwriting, gesture recognition, part-of-speech tagging, musical score following, partial discharges and bioinformatics.

Download sample and source code
Download the Accord.NET Framework

The code presented here is also part of the Accord.NET Framework. The Accord.NET Framework is a framework for developing machine learning, computer vision, computer audition, statistics and math applications. It is based on the already excellent AForge.NET Framework. Please see the starting guide for mode details. The latest version of the framework includes the latest version of this code plus many other statistics and machine learning tools.

Introduction

Hidden Markov Models were first described in a series of statistical papers by Leonard E. Baum and other authors in the second half of the 1960s. One of the first applications of HMMs was speech recognition, starting in the mid-1970s. Indeed, one of the most comprehensive explanations on the topic was published in “A Tutorial On Hidden Markov Models And Selected Applications in Speech Recognition”, by Lawrence R. Rabiner in 1989. In the second half of the 1980s, HMMs began to be applied to the analysis of biological sequences, in particular DNA. Since then, they have become ubiquitous in the field of bioinformatics.

Dynamical systems of discrete nature assumed to be governed by a Markov chain emits a sequence of observable outputs. Under the Markov assumption, it is also assumed that the latest output depends only on the current state of the system. Such states are often not known from the observer when only the output values are observable.

Hidden Markov Models attempt to model such systems and allow, among other things, (1) to infer the most likely sequence of states that produced a given output sequence, to (2) infer which will be the most likely next state (and thus predicting the next output) and (3) calculate the probability that a given sequence of outputs originated from the system (allowing the use of hidden Markov models for sequence classification).

The “hidden” in Hidden Markov Models comes from the fact that the observer does not know in which state the system may be in, but has only a probabilistic insight on where it should be.

Definition

Hidden Markov Models can be seem as finite state machines where for each sequence unit observation there is a state transition and, for each state, there is a output symbol emission.

Notation

Traditionally, HMMs have been defined by the following quintuple:

$lambda = (N, M, A, B, pi)$

where

N is the number of states for the model
M is the number of distinct observations symbols per state, i.e. the discrete alphabet size.
A is the NxN state transition probability distribution given in the form of a matrix A = {a_ij}
B is the NxM observation symbol probability distribution given in the form of a matrix B = {b_j(k)}
π is the initial state distribution vector π = {π_i}

Note that, if we opt out the structure parameters M and N we have the more often used compact notation

$lambda = (A, B, pi)$

Canonical problems

There are three canonical problems associated with hidden Markov models, which I’ll quote from Wikipedia:

Given the parameters of the model, compute the probability of a particular output sequence. This requires summation over all possible state sequences, but can be done efficiently using the Forward algorithm, which is a form of dynamic programming.
Given the parameters of the model and a particular output sequence, find the state sequence that is most likely to have generated that output sequence. This requires finding a maximum over all possible state sequences, but can similarly be solved efficiently by the Viterbi algorithm.
Given an output sequence or a set of such sequences, find the most likely set of state transition and output probabilities. In other words, derive the maximum likelihood estimate of the parameters of the HMM given a dataset of output sequences. No tractable algorithm is known for solving this problem exactly, but a local maximum likelihood can be derived efficiently using the Baum-Welch algorithm or the Baldi-Chauvin algorithm. The Baum-Welch algorithm is an example of a forward-backward algorithm, and is a special case of the Expectation-maximization algorithm.

The solution for those problems are exactly what makes Hidden Markov Models useful. The ability to learn from the data (using the solution of problem 3) and then become able to make predictions (solution to problem 2) and able to classify sequences (solution of problem 2) is nothing but applied machine learning. From this perspective, HMMs can just be seem as supervisioned sequence classifiers and sequence predictors with some other useful interesting properties.

Choosing the structure

Choosing the structure for a hidden Markov model is not always obvious. The number of states depend on the application and to what interpretation one is willing to give to the hidden states. Some domain knowledge is required to build a suitable model and also to choose the initial parameters that an HMM can take. There is also some trial and error involved, and there are sometimes complex tradeoffs that have to be made between model complexity and difficulty of learning, just as is the case with most machine learning techniques.

Additional information can be found on http://www.cse.unsw.edu.au/~waleed/phd/tr9806/node12.html.

Algorithms

The solution to the three canonical problems are the algorithms that makes HMMs useful. Each of the three problems are described in the three subsections below.

Evaluation

The first canonical problem is the evaluation of the probability of a particular output sequence. It can be efficiently computed using either the Viterbi-forward or the Forward algorithms, both of which are forms of dynamic programming.

The Viterbi algorithm originally computes the most likely sequence of states which has originated a sequence of observations. In doing so, it is also able to return the probability of traversing this particular sequence of states. So to obtain Viterbi probabilities, please refer to the Decoding problem referred below.

The Forward algorithm, unlike the Viterbi algorithm, does not find a particular sequence of states; instead it computes the probability that any sequence of states has produced the sequence of observations. In both algorithms, a matrix is used to store computations about the possible state sequence paths that the model can assume. The forward algorithm also plays a key role in the Learning problem, and is thus implemented as a separate method.

<!-- code formatted by http://manoli.net/csharpformat/ -->
<pre>

<span>/// &lt;summary&gt;</span>
<span>///   Calculates the probability that this model has generated the given sequence.</span>
<span>/// &lt;/summary&gt;</span>
<span>/// &lt;remarks&gt;</span>
<span>///   Evaluation problem. Given the HMM  M = (A, B, pi) and  the observation</span>
<span>///   sequence O = {o1, o2, ..., oK}, calculate the probability that model</span>
<span>///   M has generated sequence O. This can be computed efficiently using the</span>
<span>///   either the Viterbi or the Forward algorithms.</span>
<span>/// &lt;/remarks&gt;</span>
<span>/// &lt;param name="observations"&gt;</span>
<span>///   A sequence of observations.</span>
<span>/// &lt;/param&gt;</span>
<span>/// &lt;param name="logarithm"&gt;</span>
<span>///   True to return the log-likelihood, false to return</span>
<span>///   the likelihood. Default is false.</span>
<span>/// &lt;/param&gt;</span>
<span>/// &lt;returns&gt;</span>
<span>///   The probability that the given sequence has been generated by this model.</span>
<span>/// &lt;/returns&gt;</span>
<span>public</span> <span>double</span> Evaluate(<span>int</span>[] observations, <span>bool</span> logarithm)
{
    <span>if</span> (observations == <span>null</span>)
        <span>throw</span> <span>new</span> ArgumentNullException(<span>"observations"</span>);

    <span>if</span> (observations.Length == 0)
        <span>return</span> 0.0;


    <span>// Forward algorithm</span>
    <span>double</span> likelihood = 0;
    <span>double</span>[] coefficients;

    <span>// Compute forward probabilities</span>
    forward(observations, <span>out</span> coefficients);

    <span>for</span> (<span>int</span> i = 0; i &lt; coefficients.Length; i++)
        likelihood += Math.Log(coefficients[i]);

    <span>// Return the sequence probability</span>
    <span>return</span> logarithm ? likelihood : Math.Exp(likelihood);
}

<span>/// &lt;summary&gt;</span>
<span>///   Baum-Welch forward pass (with scaling)</span>
<span>/// &lt;/summary&gt;</span>
<span>/// &lt;remarks&gt;</span>
<span>///   Reference: http://courses.media.mit.edu/2010fall/mas622j/ProblemSets/ps4/tutorial.pdf</span>
<span>/// &lt;/remarks&gt;</span>
<span>private</span> <span>double</span>[,] forward(<span>int</span>[] observations, <span>out</span> <span>double</span>[] c)
{
    <span>int</span> T = observations.Length;
    <span>double</span>[] pi = Probabilities;
    <span>double</span>[,] A = Transitions;

    <span>double</span>[,] fwd = <span>new</span> <span>double</span>[T, States];
    c = <span>new</span> <span>double</span>[T];


    <span>// 1. Initialization</span>
    <span>for</span> (<span>int</span> i = 0; i &lt; States; i++)
        c[0] += fwd[0, i] = pi[i] * B[i, observations[0]];

    <span>if</span> (c[0] != 0) <span>// Scaling</span>
    {
        <span>for</span> (<span>int</span> i = 0; i &lt; States; i++)
            fwd[0, i] = fwd[0, i] / c[0];
    }


    <span>// 2. Induction</span>
    <span>for</span> (<span>int</span> t = 1; t &lt; T; t++)
    {
        <span>for</span> (<span>int</span> i = 0; i &lt; States; i++)
        {
            <span>double</span> p = B[i, observations[t]];

            <span>double</span> sum = 0.0;
            <span>for</span> (<span>int</span> j = 0; j &lt; States; j++)
                sum += fwd[t - 1, j] * A[j, i];
            fwd[t, i] = sum * p;

            c[t] += fwd[t, i]; <span>// scaling coefficient</span>
        }

        <span>if</span> (c[t] != 0) <span>// Scaling</span>
        {
            <span>for</span> (<span>int</span> i = 0; i &lt; States; i++)
                fwd[t, i] = fwd[t, i] / c[t];
        }
    }

    <span>return</span> fwd;
}

<pre>

/// <summary>

/// Calculates the probability that this model has generated the given sequence.

/// </summary>

/// <remarks>

/// Evaluation problem. Given the HMM M = (A, B, pi) and the observation

/// sequence O = {o1, o2, ..., oK}, calculate the probability that model

/// M has generated sequence O. This can be computed efficiently using the

/// either the Viterbi or the Forward algorithms.

/// </remarks>

/// <param name="observations">

/// A sequence of observations.

/// </param>

/// <param name="logarithm">

/// True to return the log-likelihood, false to return

/// the likelihood. Default is false.

/// </param>

/// <returns>

/// The probability that the given sequence has been generated by this model.

/// </returns>

public double Evaluate(int[] observations, bool logarithm)

{

if (observations == null)

throw new ArgumentNullException("observations");

if (observations.Length == 0)

return 0.0;

// Forward algorithm

double likelihood = 0;

double[] coefficients;

// Compute forward probabilities

forward(observations, out coefficients);

for (int i = 0; i < coefficients.Length; i++)

likelihood += Math.Log(coefficients[i]);

// Return the sequence probability

return logarithm ? likelihood : Math.Exp(likelihood);

}

/// <summary>

/// Baum-Welch forward pass (with scaling)

/// </summary>

/// <remarks>

/// Reference: http://courses.media.mit.edu/2010fall/mas622j/ProblemSets/ps4/tutorial.pdf

/// </remarks>

private double[,] forward(int[] observations, out double[] c)

{

int T = observations.Length;

double[] pi = Probabilities;

double[,] A = Transitions;

double[,] fwd = new double[T, States];

c = new double[T];

// 1. Initialization

for (int i = 0; i < States; i++)

c[0] += fwd[0, i] = pi[i] * B[i, observations[0]];

if (c[0] != 0) // Scaling

{

for (int i = 0; i < States; i++)

fwd[0, i] = fwd[0, i] / c[0];

}

// 2. Induction

for (int t = 1; t < T; t++)

{

for (int i = 0; i < States; i++)

{

double p = B[i, observations[t]];

double sum = 0.0;

for (int j = 0; j < States; j++)

sum += fwd[t - 1, j] * A[j, i];

fwd[t, i] = sum * p;

c[t] += fwd[t, i]; // scaling coefficient

}

if (c[t] != 0) // Scaling

{

for (int i = 0; i < States; i++)

fwd[t, i] = fwd[t, i] / c[t];

}

return fwd;

}

Decoding

The second canonical problem is the discovery of the most likely sequence of states that generated a given output sequence. This can be computed efficiently using the Viterbi algorithm. A trackback is used to detect the maximum probability path travelled by the algorithm. The probability of travelling such sequence is also computed in the process.

<span>/// &lt;summary&gt;</span><br /><span>///   Calculates the most likely sequence of hidden states</span><br /><span>///   that produced the given observation sequence.</span><br /><span>/// &lt;/summary&gt;</span><br /><span>/// &lt;remarks&gt;</span><br /><span>///   Decoding problem. Given the HMM M = (A, B, pi) and  the observation sequence </span><br /><span>///   O = {o1,o2, ..., oK}, calculate the most likely sequence of hidden states Si</span><br /><span>///   that produced this observation sequence O. This can be computed efficiently</span><br /><span>///   using the Viterbi algorithm.</span><br /><span>/// &lt;/remarks&gt;</span><br /><span>/// &lt;param name=&quot;observations&quot;&gt;A sequence of observations.&lt;/param&gt;</span><br /><span>/// &lt;param name=&quot;probability&quot;&gt;The state optimized probability.&lt;/param&gt;</span><br /><span>/// &lt;returns&gt;The sequence of states that most likely produced the sequence.&lt;/returns&gt;</span><br /><span>public</span> <span>int</span>[] Decode(<span>int</span>[] observations, <span>out</span> <span>double</span> probability)<br />{<br />    <span>// Viterbi algorithm.</span><br /><br />    <span>int</span> T = observations.Length;<br />    <span>int</span> states = States;<br />    <span>int</span> minState;<br />    <span>double</span> minWeight;<br />    <span>double</span> weight;<br /><br />    <span>int</span>[,] s = <span>new</span> <span>int</span>[states, T];<br />    <span>double</span>[,] a = <span>new</span> <span>double</span>[states, T];<br /><br /><br />    <span>// Base</span><br />    <span>for</span> (<span>int</span> i = 0; i &lt; states; i++)<br />    {<br />        a[i, 0] = (-1.0 * System.Math.Log(pi[i])) - System.Math.Log(B[i, observations[0]]);<br />    }<br /><br />    <span>// Induction</span><br />    <span>for</span> (<span>int</span> t = 1; t &lt; T; t++)<br />    {<br />        <span>for</span> (<span>int</span> j = 0; j &lt; states; j++)<br />        {<br />            minState = 0;<br />            minWeight = a[0, t - 1] - System.Math.Log(A[0, j]);<br /><br />            <span>for</span> (<span>int</span> i = 1; i &lt; states; i++)<br />            {<br />                weight = a[i, t - 1] - System.Math.Log(A[i, j]);<br /><br />                <span>if</span> (weight &lt; minWeight)<br />                {<br />                    minState = i;<br />                    minWeight = weight;<br />                }<br />            }<br /><br />            a[j, t] = minWeight - System.Math.Log(B[j, observations[t]]);<br />            s[j, t] = minState;<br />        }<br />    }<br /><br /><br />    <span>// Find minimum value for time T-1</span><br />    minState = 0;<br />    minWeight = a[0, T - 1];<br /><br />    <span>for</span> (<span>int</span> i = 1; i &lt; states; i++)<br />    {<br />        <span>if</span> (a[i, T - 1] &lt; minWeight)<br />        {<br />            minState = i;<br />            minWeight = a[i, T - 1];<br />        }<br />    }<br /><br />    <span>// Trackback</span><br />    <span>int</span>[] path = <span>new</span> <span>int</span>[T];<br />    path[T - 1] = minState;<br /><br />    <span>for</span> (<span>int</span> t = T - 2; t &gt;= 0; t--)<br />        path[t] = s[path[t + 1], t + 1];<br /><br /><br />    probability = System.Math.Exp(-minWeight);<br />    <span>return</span> path;<br />}

/// <summary> /// Calculates the most likely sequence of hidden states /// that produced the given observation sequence. /// </summary> /// <remarks> /// Decoding problem. Given the HMM M = (A, B, pi) and the observation sequence /// O = {o1,o2, ..., oK}, calculate the most likely sequence of hidden states Si /// that produced this observation sequence O. This can be computed efficiently /// using the Viterbi algorithm. /// </remarks> /// <param name="observations">A sequence of observations.</param> /// <param name="probability">The state optimized probability.</param> /// <returns>The sequence of states that most likely produced the sequence.</returns> public int[] Decode(int[] observations, out double probability) { // Viterbi algorithm. int T = observations.Length; int states = States; int minState; double minWeight; double weight; int[,] s = new int[states, T]; double[,] a = new double[states, T]; // Base for (int i = 0; i < states; i++) { a[i, 0] = (-1.0 * System.Math.Log(pi[i])) - System.Math.Log(B[i, observations[0]]); } // Induction for (int t = 1; t < T; t++) { for (int j = 0; j < states; j++) { minState = 0; minWeight = a[0, t - 1] - System.Math.Log(A[0, j]); for (int i = 1; i < states; i++) { weight = a[i, t - 1] - System.Math.Log(A[i, j]); if (weight < minWeight) { minState = i; minWeight = weight; } } a[j, t] = minWeight - System.Math.Log(B[j, observations[t]]); s[j, t] = minState; } } // Find minimum value for time T-1 minState = 0; minWeight = a[0, T - 1]; for (int i = 1; i < states; i++) { if (a[i, T - 1] < minWeight) { minState = i; minWeight = a[i, T - 1]; } } // Trackback int[] path = new int[T]; path[T - 1] = minState; for (int t = T - 2; t >= 0; t--) path[t] = s[path[t + 1], t + 1]; probability = System.Math.Exp(-minWeight); return path; }

Learning

The third and last problem is the problem of learning the most likely parameters that best models a system given a set of sequences originated from this system. Most implementations I’ve seem did not consider the problem of learning from a set of sequences, but only from a single sequence at a time. The algorithm below, however, is fully suitable to learn from a set of sequences and also uses scaling, which is another thing I have not seem in other implementations.

The source code follows the original algorithm by Rabiner (1989). There are, however, some known issues with the algorithms detailed in Rabiner’s paper. More information about those issues is available in a next section of this article entitled “Remarks”.

<!-- code formatted by http://manoli.net/csharpformat/ -->
<pre>
<span>/// &lt;summary&gt;</span>
<span>///   Runs the Baum-Welch learning algorithm for hidden Markov models.</span>
<span>/// &lt;/summary&gt;</span>
<span>/// &lt;remarks&gt;</span>
<span>///   Learning problem. Given some training observation sequences O = {o1, o2, ..., oK}</span>
<span>///   and general structure of HMM (numbers of hidden and visible states), determine</span>
<span>///   HMM parameters M = (A, B, pi) that best fit training data. </span>
<span>/// &lt;/remarks&gt;</span>
<span>/// &lt;param name="iterations"&gt;</span>
<span>///   The maximum number of iterations to be performed by the learning algorithm. If</span>
<span>///   specified as zero, the algorithm will learn until convergence of the model average</span>
<span>///   likelihood respecting the desired limit.</span>
<span>/// &lt;/param&gt;</span>
<span>/// &lt;param name="observations"&gt;</span>
<span>///   An array of observation sequences to be used to train the model.</span>
<span>/// &lt;/param&gt;</span>
<span>/// &lt;param name="tolerance"&gt;</span>
<span>///   The likelihood convergence limit L between two iterations of the algorithm. The</span>
<span>///   algorithm will stop when the change in the likelihood for two consecutive iterations</span>
<span>///   has not changed by more than L percent of the likelihood. If left as zero, the</span>
<span>///   algorithm will ignore this parameter and iterates over a number of fixed iterations</span>
<span>///   specified by the previous parameter.</span>
<span>/// &lt;/param&gt;</span>
<span>/// &lt;returns&gt;</span>
<span>///   The average log-likelihood for the observations after the model has been trained.</span>
<span>/// &lt;/returns&gt;</span>
<span>public</span> <span>double</span> Learn(<span>int</span>[][] observations, <span>int</span> iterations, <span>double</span> tolerance)
{
    <span>if</span> (iterations == 0 &amp;&amp; tolerance == 0)
        <span>throw</span> <span>new</span> ArgumentException(<span>"Iterations and limit cannot be both zero."</span>);

    <span>// Baum-Welch algorithm.</span>

    <span>// The Baum–Welch algorithm is a particular case of a generalized expectation-maximization</span>
    <span>// (GEM) algorithm. It can compute maximum likelihood estimates and posterior mode estimates</span>
    <span>// for the parameters (transition and emission probabilities) of an HMM, when given only</span>
    <span>// emissions as training data.</span>

    <span>// The algorithm has two steps:</span>
    <span>//  - Calculating the forward probability and the backward probability for each HMM state;</span>
    <span>//  - On the basis of this, determining the frequency of the transition-emission pair values</span>
    <span>//    and dividing it by the probability of the entire string. This amounts to calculating</span>
    <span>//    the expected count of the particular transition-emission pair. Each time a particular</span>
    <span>//    transition is found, the value of the quotient of the transition divided by the probability</span>
    <span>//    of the entire string goes up, and this value can then be made the new value of the transition.</span>


    <span>int</span> N = observations.Length;
    <span>int</span> currentIteration = 1;
    <span>bool</span> stop = <span>false</span>;

    <span>double</span>[] pi = Probabilities;
    <span>double</span>[,] A = Transitions;


    <span>// Initialization</span>
    <span>double</span>[][, ,] epsilon = <span>new</span> <span>double</span>[N][, ,]; <span>// also referred as ksi or psi</span>
    <span>double</span>[][,] gamma = <span>new</span> <span>double</span>[N][,];

    <span>for</span> (<span>int</span> i = 0; i &lt; N; i++)
    {
        <span>int</span> T = observations[i].Length;
        epsilon[i] = <span>new</span> <span>double</span>[T, States, States];
        gamma[i] = <span>new</span> <span>double</span>[T, States];
    }


    <span>// Calculate initial model log-likelihood</span>
    <span>double</span> oldLikelihood = Double.MinValue;
    <span>double</span> newLikelihood = 0;


    <span>do</span> <span>// Until convergence or max iterations is reached</span>
    {
        <span>// For each sequence in the observations input</span>
        <span>for</span> (<span>int</span> i = 0; i &lt; N; i++)
        {
            var sequence = observations[i];
            <span>int</span> T = sequence.Length;
            <span>double</span>[] scaling;

            <span>// 1st step - Calculating the forward probability and the</span>
            <span>//            backward probability for each HMM state.</span>
            <span>double</span>[,] fwd = forward(observations[i], <span>out</span> scaling);
            <span>double</span>[,] bwd = backward(observations[i], scaling);


            <span>// 2nd step - Determining the frequency of the transition-emission pair values</span>
            <span>//            and dividing it by the probability of the entire string.</span>


            <span>// Calculate gamma values for next computations</span>
            <span>for</span> (<span>int</span> t = 0; t &lt; T; t++)
            {
                <span>double</span> s = 0;

                <span>for</span> (<span>int</span> k = 0; k &lt; States; k++)
                    s += gamma[i][t, k] = fwd[t, k] * bwd[t, k];

                <span>if</span> (s != 0) <span>// Scaling</span>
                {
                    <span>for</span> (<span>int</span> k = 0; k &lt; States; k++)
                        gamma[i][t, k] /= s;
                }
            }

            <span>// Calculate epsilon values for next computations</span>
            <span>for</span> (<span>int</span> t = 0; t &lt; T - 1; t++)
            {
                <span>double</span> s = 0;

                <span>for</span> (<span>int</span> k = 0; k &lt; States; k++)
                    <span>for</span> (<span>int</span> l = 0; l &lt; States; l++)
                        s += epsilon[i][t, k, l] = fwd[t, k] * A[k, l] * bwd[t + 1, l] * B[l, sequence[t + 1]];

                <span>if</span> (s != 0) <span>// Scaling</span>
                {
                    <span>for</span> (<span>int</span> k = 0; k &lt; States; k++)
                        <span>for</span> (<span>int</span> l = 0; l &lt; States; l++)
                            epsilon[i][t, k, l] /= s;
                }
            }

            <span>// Compute log-likelihood for the given sequence</span>
            <span>for</span> (<span>int</span> t = 0; t &lt; scaling.Length; t++)
                newLikelihood += Math.Log(scaling[t]);
        }


        <span>// Average the likelihood for all sequences</span>
        newLikelihood /= observations.Length;


        <span>// Check if the model has converged or we should stop</span>
        <span>if</span> (checkConvergence(oldLikelihood, newLikelihood,
            currentIteration, iterations, tolerance))
        {
            stop = <span>true</span>;
        }

        <span>else</span>
        {
            <span>// 3. Continue with parameter re-estimation</span>
            currentIteration++;
            oldLikelihood = newLikelihood;
            newLikelihood = 0.0;


            <span>// 3.1 Re-estimation of initial state probabilities </span>
            <span>for</span> (<span>int</span> k = 0; k &lt; States; k++)
            {
                <span>double</span> sum = 0;
                <span>for</span> (<span>int</span> i = 0; i &lt; N; i++)
                    sum += gamma[i][0, k];
                pi[k] = sum / N;
            }

            <span>// 3.2 Re-estimation of transition probabilities </span>
            <span>for</span> (<span>int</span> i = 0; i &lt; States; i++)
            {
                <span>for</span> (<span>int</span> j = 0; j &lt; States; j++)
                {
                    <span>double</span> den = 0, num = 0;

                    <span>for</span> (<span>int</span> k = 0; k &lt; N; k++)
                    {
                        <span>int</span> T = observations[k].Length;

                        <span>for</span> (<span>int</span> l = 0; l &lt; T - 1; l++)
                            num += epsilon[k][l, i, j];

                        <span>for</span> (<span>int</span> l = 0; l &lt; T - 1; l++)
                            den += gamma[k][l, i];
                    }

                    A[i, j] = (den != 0) ? num / den : 0.0;
                }
            }

            <span>// 3.3 Re-estimation of emission probabilities</span>
            <span>for</span> (<span>int</span> i = 0; i &lt; States; i++)
            {
                <span>for</span> (<span>int</span> j = 0; j &lt; Symbols; j++)
                {
                    <span>double</span> den = 0, num = 0;

                    <span>for</span> (<span>int</span> k = 0; k &lt; N; k++)
                    {
                        <span>int</span> T = observations[k].Length;

                        <span>for</span> (<span>int</span> l = 0; l &lt; T; l++)
                        {
                            <span>if</span> (observations[k][l] == j)
                                num += gamma[k][l, i];
                        }

                        <span>for</span> (<span>int</span> l = 0; l &lt; T; l++)
                            den += gamma[k][l, i];
                    }

                    <span>// avoid locking a parameter in zero.</span>
                    B[i, j] = (num == 0) ? 1e-10 : num / den;
                }
            }

        }

    } <span>while</span> (!stop);


    <span>// Returns the model average log-likelihood</span>
    <span>return</span> newLikelihood;
}

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

<pre>

/// <summary>

/// Runs the Baum-Welch learning algorithm for hidden Markov models.

/// </summary>

/// <remarks>

/// Learning problem. Given some training observation sequences O = {o1, o2, ..., oK}

/// and general structure of HMM (numbers of hidden and visible states), determine

/// HMM parameters M = (A, B, pi) that best fit training data.

/// </remarks>

/// <param name="iterations">

/// The maximum number of iterations to be performed by the learning algorithm. If

/// specified as zero, the algorithm will learn until convergence of the model average

/// likelihood respecting the desired limit.

/// </param>

/// <param name="observations">

/// An array of observation sequences to be used to train the model.

/// </param>

/// <param name="tolerance">

/// The likelihood convergence limit L between two iterations of the algorithm. The

/// algorithm will stop when the change in the likelihood for two consecutive iterations

/// has not changed by more than L percent of the likelihood. If left as zero, the

/// algorithm will ignore this parameter and iterates over a number of fixed iterations

/// specified by the previous parameter.

/// </param>

/// <returns>

/// The average log-likelihood for the observations after the model has been trained.

/// </returns>

public double Learn(int[][] observations, int iterations, double tolerance)

{

if (iterations == 0 && tolerance == 0)

throw new ArgumentException("Iterations and limit cannot be both zero.");

// Baum-Welch algorithm.

// The Baum–Welch algorithm is a particular case of a generalized expectation-maximization

// (GEM) algorithm. It can compute maximum likelihood estimates and posterior mode estimates

// for the parameters (transition and emission probabilities) of an HMM, when given only

// emissions as training data.

// The algorithm has two steps:

// - Calculating the forward probability and the backward probability for each HMM state;

// - On the basis of this, determining the frequency of the transition-emission pair values

// and dividing it by the probability of the entire string. This amounts to calculating

// the expected count of the particular transition-emission pair. Each time a particular

// transition is found, the value of the quotient of the transition divided by the probability

// of the entire string goes up, and this value can then be made the new value of the transition.

int N = observations.Length;

int currentIteration = 1;

bool stop = false;

double[] pi = Probabilities;

double[,] A = Transitions;

// Initialization

double[][, ,] epsilon = new double[N][, ,]; // also referred as ksi or psi

double[][,] gamma = new double[N][,];

for (int i = 0; i < N; i++)

{

int T = observations[i].Length;

epsilon[i] = new double[T, States, States];

gamma[i] = new double[T, States];

}

// Calculate initial model log-likelihood

double oldLikelihood = Double.MinValue;

double newLikelihood = 0;

do // Until convergence or max iterations is reached

{

// For each sequence in the observations input

for (int i = 0; i < N; i++)

{

var sequence = observations[i];

int T = sequence.Length;

double[] scaling;

// 1st step - Calculating the forward probability and the

// backward probability for each HMM state.

double[,] fwd = forward(observations[i], out scaling);

double[,] bwd = backward(observations[i], scaling);

// 2nd step - Determining the frequency of the transition-emission pair values

// and dividing it by the probability of the entire string.

// Calculate gamma values for next computations

for (int t = 0; t < T; t++)

{

double s = 0;

for (int k = 0; k < States; k++)

s += gamma[i][t, k] = fwd[t, k] * bwd[t, k];

if (s != 0) // Scaling

{

for (int k = 0; k < States; k++)

gamma[i][t, k] /= s;

}

// Calculate epsilon values for next computations

for (int t = 0; t < T - 1; t++)

{

double s = 0;

for (int k = 0; k < States; k++)

for (int l = 0; l < States; l++)

s += epsilon[i][t, k, l] = fwd[t, k] * A[k, l] * bwd[t + 1, l] * B[l, sequence[t + 1]];

if (s != 0) // Scaling

{

for (int k = 0; k < States; k++)

for (int l = 0; l < States; l++)

epsilon[i][t, k, l] /= s;

}

// Compute log-likelihood for the given sequence

for (int t = 0; t < scaling.Length; t++)

newLikelihood += Math.Log(scaling[t]);

}

// Average the likelihood for all sequences

newLikelihood /= observations.Length;

// Check if the model has converged or we should stop

if (checkConvergence(oldLikelihood, newLikelihood,

currentIteration, iterations, tolerance))

{

stop = true;

}

else

{

// 3. Continue with parameter re-estimation

currentIteration++;

oldLikelihood = newLikelihood;

newLikelihood = 0.0;

// 3.1 Re-estimation of initial state probabilities

for (int k = 0; k < States; k++)

{

double sum = 0;

for (int i = 0; i < N; i++)

sum += gamma[i][0, k];

pi[k] = sum / N;

}

// 3.2 Re-estimation of transition probabilities

for (int i = 0; i < States; i++)

{

for (int j = 0; j < States; j++)

{

double den = 0, num = 0;

for (int k = 0; k < N; k++)

{

int T = observations[k].Length;

for (int l = 0; l < T - 1; l++)

num += epsilon[k][l, i, j];

for (int l = 0; l < T - 1; l++)

den += gamma[k][l, i];

}

A[i, j] = (den != 0) ? num / den : 0.0;

}

// 3.3 Re-estimation of emission probabilities

for (int i = 0; i < States; i++)

{

for (int j = 0; j < Symbols; j++)

{

double den = 0, num = 0;

for (int k = 0; k < N; k++)

{

int T = observations[k].Length;

for (int l = 0; l < T; l++)

{

if (observations[k][l] == j)

num += gamma[k][l, i];

}

for (int l = 0; l < T; l++)

den += gamma[k][l, i];

}

// avoid locking a parameter in zero.

B[i, j] = (num == 0) ? 1e-10 : num / den;

}

} while (!stop);

// Returns the model average log-likelihood

return newLikelihood;

}

Using the code

Lets suppose we have gathered some sequences from a system we wish to model. The sequences are expressed as a integer array such as:

<span>int</span>[][] sequences = <span>new</span> <span>int</span>[][] <br />{<br />    <span>new</span> <span>int</span>[] { 0,1,1,1,1,1,1 },<br />    <span>new</span> <span>int</span>[] { 0,1,1,1 },<br />    <span>new</span> <span>int</span>[] { 0,1,1,1,1 },<br />    <span>new</span> <span>int</span>[] { 0,1, },<br />    <span>new</span> <span>int</span>[] { 0,1,1 },<br />};

int[][] sequences = new int[][] { new int[] { 0,1,1,1,1,1,1 }, new int[] { 0,1,1,1 }, new int[] { 0,1,1,1,1 }, new int[] { 0,1, }, new int[] { 0,1,1 }, };

For us, it can be obvious to see that the system is outputting sequences that always start with a zero and have one or more ones at the end. But lets try to fit a Hidden Markov Model to predict those sequences.

<span>// Creates a new Hidden Markov Model with 2 states for</span><br /><span>//  an output alphabet of two characters (zero and one)</span><br />HiddenMarkovModel hmm = <span>new</span> HiddenMarkovModel(2, 2);<br /><br /><span>// Try to fit the model to the data until the difference in</span><br /><span>//  the average likelihood changes only by as little as 0.01</span><br />hmm.Learn(sequences, 0.01);

// Creates a new Hidden Markov Model with 2 states for // an output alphabet of two characters (zero and one) HiddenMarkovModel hmm = new HiddenMarkovModel(2, 2); // Try to fit the model to the data until the difference in // the average likelihood changes only by as little as 0.01 hmm.Learn(sequences, 0.01);

Once the model is trained, lets test to see if it recognizes some sequences:

<span>// Calculate the probability that the given</span><br /><span>//  sequences originated from the model</span><br /><span>double</span> l1 = hmm.Evaluate(<span>new</span> <span>int</span>[] { 0, 1 });       <span>// l1 = 0.9999</span><br /><span>double</span> l2 = hmm.Evaluate(<span>new</span> <span>int</span>[] { 0, 1, 1, 1 }); <span>// l2 = 0.9999</span><br /><br /><span>double</span> l3 = hmm.Evaluate(<span>new</span> <span>int</span>[] { 1, 1 });       <span>// l3 = 0.0000</span><br /><span>double</span> l4 = hmm.Evaluate(<span>new</span> <span>int</span>[] { 1, 0, 0, 0 }); <span>// l4 = 0.0000</span>

// Calculate the probability that the given // sequences originated from the model double l1 = hmm.Evaluate(new int[] { 0, 1 }); // l1 = 0.9999 double l2 = hmm.Evaluate(new int[] { 0, 1, 1, 1 }); // l2 = 0.9999 double l3 = hmm.Evaluate(new int[] { 1, 1 }); // l3 = 0.0000 double l4 = hmm.Evaluate(new int[] { 1, 0, 0, 0 }); // l4 = 0.0000

Of course the model performs well as this a rather simple example. A more useful test case would consist of allowing for some errors in the input sequences in the hope that the model will become more tolerant to measurement errors.

<span>int</span>[][] sequences = <span>new</span> <span>int</span>[][] <br />{<br />    <span>new</span> <span>int</span>[] { 0,1,1,1,1,0,1,1,1,1 },<br />    <span>new</span> <span>int</span>[] { 0,1,1,1,0,1,1,1,1,1 },<br />    <span>new</span> <span>int</span>[] { 0,1,1,1,1,1,1,1,1,1 },<br />    <span>new</span> <span>int</span>[] { 0,1,1,1,1,1         },<br />    <span>new</span> <span>int</span>[] { 0,1,1,1,1,1,1       },<br />    <span>new</span> <span>int</span>[] { 0,1,1,1,1,1,1,1,1,1 },<br />    <span>new</span> <span>int</span>[] { 0,1,1,1,1,1,1,1,1,1 },<br />};<br /><br /><span>// Creates a new Hidden Markov Model with 3 states for</span><br /><span>//  an output alphabet of two characters (zero and one)</span><br />HiddenMarkovModel hmm = <span>new</span> HiddenMarkovModel(2, 3);<br /><br /><span>// Try to fit the model to the data until the difference in</span><br /><span>//  the average likelihood changes only by as little as 0.0001</span><br />hmm.Learn(sequences, 0.0001);<br /><br /><span>// Calculate the probability that the given</span><br /><span>//  sequences originated from the model</span><br /><span>double</span> l1 = hmm.Evaluate(<span>new</span> <span>int</span>[] { 0,1 });      <span>// 0.9999  </span><br /><span>double</span> l2 = hmm.Evaluate(<span>new</span> <span>int</span>[] { 0,1,1,1 });  <span>// 0.9166</span><br /><br /><span>double</span> l3 = hmm.Evaluate(<span>new</span> <span>int</span>[] { 1,1 });      <span>// 0.0000</span><br /><span>double</span> l4 = hmm.Evaluate(<span>new</span> <span>int</span>[] { 1,0,0,0 });  <span>// 0.0000</span><br /><br /><span>double</span> l5 = hmm.Evaluate(<span>new</span> <span>int</span>[] { 0,1,0,1,1,1,1,1,1 }); <span>// 0.0342</span><br /><span>double</span> l6 = hmm.Evaluate(<span>new</span> <span>int</span>[] { 0,1,1,1,1,1,1,0,1 }); <span>// 0.0342</span>

int[][] sequences = new int[][] { new int[] { 0,1,1,1,1,0,1,1,1,1 }, new int[] { 0,1,1,1,0,1,1,1,1,1 }, new int[] { 0,1,1,1,1,1,1,1,1,1 }, new int[] { 0,1,1,1,1,1 }, new int[] { 0,1,1,1,1,1,1 }, new int[] { 0,1,1,1,1,1,1,1,1,1 }, new int[] { 0,1,1,1,1,1,1,1,1,1 }, }; // Creates a new Hidden Markov Model with 3 states for // an output alphabet of two characters (zero and one) HiddenMarkovModel hmm = new HiddenMarkovModel(2, 3); // Try to fit the model to the data until the difference in // the average likelihood changes only by as little as 0.0001 hmm.Learn(sequences, 0.0001); // Calculate the probability that the given // sequences originated from the model double l1 = hmm.Evaluate(new int[] { 0,1 }); // 0.9999 double l2 = hmm.Evaluate(new int[] { 0,1,1,1 }); // 0.9166 double l3 = hmm.Evaluate(new int[] { 1,1 }); // 0.0000 double l4 = hmm.Evaluate(new int[] { 1,0,0,0 }); // 0.0000 double l5 = hmm.Evaluate(new int[] { 0,1,0,1,1,1,1,1,1 }); // 0.0342 double l6 = hmm.Evaluate(new int[] { 0,1,1,1,1,1,1,0,1 }); // 0.0342

We can see that, despite having a very low probability, the likelihood values for the sequences containing a simulated measurement error are greater than the likelihoods for the sequences which do not follow the sequence structure at all.

In a subsequent article, we will see that those low values for the likelihoods will not be a problem because HMMs are often used in sets to form sequence classifiers. When used in such configurations, what really matters is which HMM returns the highest probability among others in the set.

Remarks

A practical issue in the use of Hidden Markov Models to model long sequences is the numerical scaling of conditional probabilities. The probability of observing a long sequence given most models is extremely small, and the use of these extremely small numbers in computations often leads to numerical instability, making application of HMMs to genome length sequences quite challenging.

There are two common approaches to dealing with small conditional probabilities. One approach is to rescale the conditional probabilities using carefully designed scaling factors, and the other approach is to work with the logarithms of the conditional probabilities. For more information on using logarithms please see the work entitled “Numerically Stable Hidden Markov Model Implementation”, by Tobias P. Mann.

Known issues

The code on this article is based on the Tutorial by Rabiner. There are, however, some problems with the scaling and other algorithms. An errata depicting all issues is available in the website “An Erratum for ‘A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition’” and is maintained by Ali Rahimi. I have not yet verified if the implementation presented here also suffers from the same mistakes explained there. This code has worked well under many situations, but I cannot guarantee its perfectness. Please use at your own risk.

Acknowledgements

Thanks to Guilherme C. Pedroso, for the help with the Baum-Welch generalization for multiple input sequences. He has also co-written a very interesting article using hidden Markov models for gesture recognition, entitled “Automatic Recognition of Finger Spelling for LIBRAS based on a Two-Layer Architecture” published in the 25th Symposium On Applied Computing (ACM SAC 2010).

References

L. R. Rabiner, "A Tutorial on Hidden Markov Models, and Selected Applications in Speech Recognition," Proc. IEEE, Vol. 77, No. 2, pp. 257–286, Feb. 1989.
Warakagoda, Narada D., “A Hybrid ANN-HMM ASR system with NN based adaptive preprocessing”, Norges Tekniske Høgskole, 2006.
A. Rahimim, “An erratum for: A tutorial on hidden Markov models and selected applications in speech recognition”, World Wide Web. Available in: http://xenia.media.mit.edu/rahimi/rabiner/rabiner-errata/rabiner-errata.html.
T. P. Mann, "Numerically stable hidden markov model implementation," World Wide Web. Available in: http://www.phrap.org/compbio/mbt599/hmm_scaling_revised.pdf
Wikipedia contributors, “Hidden Markov model,” Wikipedia, The Free Encyclopedia, http://en.wikipedia.org/wiki/Hidden_Markov_model (accessed March, 2010).

CodeProject

Kernel Functions for Machine Learning Applications

March 17, 2010 / cesarsouza / 63 Comments

In recent years, Kernel methods have received major attention, particularly due to the increased popularity of the Support Vector Machines. Kernel functions can be used in many applications as they provide a simple bridge from linearity to non-linearity for algorithms which can be expressed in terms of dot products. In this article, we will list a few kernel functions and some of their properties.

Check the source code for all kernel functions here.

Many of these functions have been incorporated in Accord.NET, a framework for creating machine learning, statistics, and computer vision applications.

Kernel Methods
Kernel Functions
Source code
See also
References

Kernel Methods

Kernel methods are a class of algorithms for pattern analysis or recognition, whose best known element is the support vector machine (SVM). The general task of pattern analysis is to find and study general types of relations (such as clusters, rankings, principal components, correlations, classifications) in general types of data (such as sequences, text documents, sets of points, vectors, images, graphs, etc) (Wikipedia, 2010a).

The main characteristic of Kernel Methods, however, is their distinct approach to this problem. Kernel methods map the data into higher dimensional spaces in the hope that in this higher-dimensional space the data could become more easily separated or better structured. There are also no constraints on the form of this mapping, which could even lead to infinite-dimensional spaces. This mapping function, however, hardly needs to be computed because of a tool called the kernel trick.

The Kernel trick

The Kernel trick is a very interesting and powerful tool. It is powerful because it provides a bridge from linearity to non-linearity to any algorithm that can expressed solely on terms of dot products between two vectors. It comes from the fact that, if we first map our input data into a higher-dimensional space, a linear algorithm operating in this space will behave non-linearly in the original input space.

K(x,y) = <φ(x),φ(y)>

Kernel Properties

Kernel functions must be continuous, symmetric, and most preferably should have a positive (semi-) definite Gram matrix. Kernels which are said to satisfy the Mercer’s theorem are positive semi-definite, meaning their kernel matrices have only non-negative Eigen values. The use of a positive definite kernel insures that the optimization problem will be convex and solution will be unique.

However, many kernel functions which aren’t strictly positive definite also have been shown to perform very well in practice. An example is the Sigmoid kernel, which, despite its wide use, it is not positive semi-definite for certain values of its parameters. Boughorbel (2005) also experimentally demonstrated that Kernels which are only conditionally positive definite can possibly outperform most classical kernels in some applications.

Kernels also can be classified as anisotropic stationary, isotropic stationary, compactly supported, locally stationary, nonstationary or separable nonstationary. Moreover, kernels can also be labeled scale-invariant or scale-dependant, which is an interesting property as scale-invariant kernels drive the training process invariant to a scaling of the data.

Choosing the Right Kernel

Choosing the most appropriate kernel highly depends on the problem at hand – and fine tuning its parameters can easily become a tedious and cumbersome task. Automatic kernel selection is possible and is discussed in the works by Tom Howley and Michael Madden.

The choice of a Kernel depends on the problem at hand because it depends on what we are trying to model. A polynomial kernel, for example, allows us to model feature conjunctions up to the order of the polynomial. Radial basis functions allows to pick out circles (or hyperspheres) – in constrast with the Linear kernel, which allows only to pick out lines (or hyperplanes).

The motivation behind the choice of a particular kernel can be very intuitive and straightforward depending on what kind of information we are expecting to extract about the data. Please see the final notes on this topic from Introduction to Information Retrieval, by Manning, Raghavan and Schütze for a better explanation on the subject.

Kernel Functions

Below is a list of some kernel functions available from the existing literature. As was the case with previous articles, every LaTeX notation for the formulas below are readily available from their alternate text html tag. I can not guarantee all of them are perfectly correct, thus use them at your own risk. Most of them have links to articles where they have been originally used or proposed.

1. Linear Kernel

Linear kernel documentation – linear kernel source code – how to create SVMs in .NET with Accord.NET

The Linear kernel is the simplest kernel function. It is given by the inner product <x,y> plus an optional constant c. Kernel algorithms using a linear kernel are often equivalent to their non-kernel counterparts, i.e. KPCA with linear kernel is the same as standard PCA.

$k(x, y) = x^T y + c$

2. Polynomial Kernel

The Polynomial kernel is a non-stationary kernel. Polynomial kernels are well suited for problems where all the training data is normalized.

$k(x, y) = (alpha x^T y + c)^d$
Adjustable parameters are the slope alpha, the constant term c and the polynomial degree d.

3. Gaussian Kernel

The Gaussian kernel is an example of radial basis function kernel.

$k(x, y) = expleft(-frac{ lVert x-y rVert ^2}{2sigma^2}right)$

Alternatively, it could also be implemented using

$k(x, y) = expleft(- gamma lVert x-y rVert ^2 )$

The adjustable parameter sigma plays a major role in the performance of the kernel, and should be carefully tuned to the problem at hand. If overestimated, the exponential will behave almost linearly and the higher-dimensional projection will start to lose its non-linear power. In the other hand, if underestimated, the function will lack regularization and the decision boundary will be highly sensitive to noise in training data.

4. Exponential Kernel

The exponential kernel is closely related to the Gaussian kernel, with only the square of the norm left out. It is also a radial basis function kernel.

$k(x, y) = expleft(-frac{ lVert x-y rVert }{2sigma^2}right)$

5. Laplacian Kernel

The Laplace Kernel is completely equivalent to the exponential kernel, except for being less sensitive for changes in the sigma parameter. Being equivalent, it is also a radial basis function kernel.

$k(x, y) = expleft(- frac{lVert x-y rVert }{sigma}right)$

It is important to note that the observations made about the sigma parameter for the Gaussian kernel also apply to the Exponential and Laplacian kernels.

6. ANOVA Kernel

The ANOVA kernel is also a radial basis function kernel, just as the Gaussian and Laplacian kernels. It is said to perform well in multidimensional regression problems (Hofmann, 2008).

$k(x, y) = sum_{k=1}^n exp (-sigma (x^k - y^k)^2)^d$

7. Hyperbolic Tangent (Sigmoid) Kernel

The Hyperbolic Tangent Kernel is also known as the Sigmoid Kernel and as the Multilayer Perceptron (MLP) kernel. The Sigmoid Kernel comes from the Neural Networks field, where the bipolar sigmoid function is often used as an activation function for artificial neurons.

$k(x, y) = tanh (alpha x^T y + c)$

It is interesting to note that a SVM model using a sigmoid kernel function is equivalent to a two-layer, perceptron neural network. This kernel was quite popular for support vector machines due to its origin from neural network theory. Also, despite being only conditionally positive definite, it has been found to perform well in practice.

There are two adjustable parameters in the sigmoid kernel, the slope alpha and the intercept constant c. A common value for alpha is 1/N, where N is the data dimension. A more detailed study on sigmoid kernels can be found in the works by Hsuan-Tien and Chih-Jen.

8. Rational Quadratic Kernel

The Rational Quadratic kernel is less computationally intensive than the Gaussian kernel and can be used as an alternative when using the Gaussian becomes too expensive.

$k(x, y) = 1 - frac{lVert x-y rVert^2}{lVert x-y rVert^2 + c}$

9. Multiquadric Kernel

The Multiquadric kernel can be used in the same situations as the Rational Quadratic kernel. As is the case with the Sigmoid kernel, it is also an example of an non-positive definite kernel.

$k(x, y) = sqrt{lVert x-y rVert^2 + c^2}$

10. Inverse Multiquadric Kernel

The Inverse Multi Quadric kernel. As with the Gaussian kernel, it results in a kernel matrix with full rank (Micchelli, 1986) and thus forms a infinite dimension feature space.

$k(x, y) = frac{1}{sqrt{lVert x-y rVert^2 + theta^2}}$

11. Circular Kernel

The circular kernel is used in geostatic applications. It is an example of an isotropic stationary kernel and is positive definite in R².

$k(x, y) = frac{2}{pi} arccos ( - frac{ lVert x-y rVert}{sigma}) - frac{2}{pi} frac{ lVert x-y rVert}{sigma} sqrt{1 - left(frac{ lVert x-y rVert}{sigma} right)^2}$
$mbox{if}~ lVert x-y rVert < sigma mbox{, zero otherwise}$

12. Spherical Kernel

The spherical kernel is similar to the circular kernel, but is positive definite in R³.

$k(x, y) = 1 - frac{3}{2} frac{lVert x-y rVert}{sigma} + frac{1}{2} left( frac{ lVert x-y rVert}{sigma} right)^3$

$mbox{if}~ lVert x-y rVert < sigma mbox{, zero otherwise}$

13. Wave Kernel

The Wave kernel is also symmetric positive semi-definite (Huang, 2008).

$k(x, y) = frac{theta}{lVert x-y rVert right} sin frac{lVert x-y rVert }{theta}$

14. Power Kernel

The Power kernel is also known as the (unrectified) triangular kernel. It is an example of scale-invariant kernel (Sahbi and Fleuret, 2004) and is also only conditionally positive definite.

$k(x,y) = - lVert x-y rVert ^d$

15. Log Kernel

The Log kernel seems to be particularly interesting for images, but is only conditionally positive definite.

$k(x,y) = - log (lVert x-y rVert ^d + 1)$

16. Spline Kernel

The Spline kernel is given as a piece-wise cubic polynomial, as derived in the works by Gunn (1998).

$k(x, y) = 1 + xy + xy~min(x,y) - frac{x+y}{2}~min(x,y)^2+frac{1}{3}min(x,y)^3$

However, what it actually mean is:

$k(x,y) = prod_{i=1}^d 1 + x_i y_i + x_i y_i min(x_i, y_i) - frac{x_i + y_i}{2} min(x_i,y_i)^2 + frac{min(x_i,y_i)^3}{3}$

With $x,y in R^d$

17. B-Spline (Radial Basis Function) Kernel

The B-Spline kernel is defined on the interval [−1, 1]. It is given by the recursive formula:

$k(x,y) = B_{2p+1}(x-y)$

$mbox{where~} p in N mbox{~with~} B_{i+1} := B_i otimes B_0.$

In the work by Bart Hamers it is given by:

$k(x, y) = prod_{p=1}^d B_{2n+1}(x_p - y_p)$

Alternatively, B_n can be computed using the explicit expression (Fomel, 2000):

$B_n(x) = frac{1}{n!} sum_{k=0}^{n+1} binom{n+1}{k} (-1)^k (x + frac{n+1}{2} - k)^n_+$

Where x₊ is defined as the truncated power function:

$x^d_+ = begin{cases} x^d, & mbox{if }x > 0 0, & mbox{otherwise} end{cases}$

18. Bessel Kernel

The Bessel kernel is well known in the theory of function spaces of fractional smoothness. It is given by:

$k(x, y) = frac{J_{v+1}( sigma lVert x-y rVert)}{ lVert x-y rVert ^ {-n(v+1)} }$

where J is the Bessel function of first kind. However, in the Kernlab for R documentation, the Bessel kernel is said to be:

$k(x,x') = - Bessel_{(nu+1)}^n (sigma |x - x'|^2)$

19. Cauchy Kernel

The Cauchy kernel comes from the Cauchy distribution (Basak, 2008). It is a long-tailed kernel and can be used to give long-range influence and sensitivity over the high dimension space.

Foreword

Contents

Introduction

Support Vector Machines

Kernel Support Vector Machines

The Kernel trick

Standard Kernels

Learning Algorithms

Sequential Minimal Optimization

Sequential Minimal Optimization for Regression

Source Code

Support Vector Machine

Kernel Support Vector Machine

Sequential Minimal Optimization

Using the code

Sample application

Classification

Regression

See also

References

Contents

Introduction

Overview

Multivariate Linear Regression in Latent Space

Algorithm 1: NIPALS

Algorithm 2: SIMPLS

Source Code

Class Diagram

Performing PLS using NIPALS

Performing PLS using SIMPLS

Multivariate Linear Regression

Using the code

Sample application

References

See also

Introduction

Hidden Markov Models

Hidden Markov Model -Based Sequence Classifier

Source Code

Computing the most likely class for a given sequence

Training each model to recognize each of the output classes

Using the Code

Sample Application

See also

References

Contents

Contents

Kernel Methods

The Kernel trick

Kernel Properties

Choosing the Right Kernel

Kernel Functions

1. Linear Kernel

2. Polynomial Kernel

3. Gaussian Kernel

4. Exponential Kernel

5. Laplacian Kernel

6. ANOVA Kernel

7. Hyperbolic Tangent (Sigmoid) Kernel

8. Rational Quadratic Kernel

9. Multiquadric Kernel

10. Inverse Multiquadric Kernel

11. Circular Kernel

12. Spherical Kernel

13. Wave Kernel

14. Power Kernel

15. Log Kernel

16. Spline Kernel

17. B-Spline (Radial Basis Function) Kernel

18. Bessel Kernel

19. Cauchy Kernel