After a preliminary article on hidden Markov models, some months ago I had finally posted the article on Hidden Conditional Random Fields (HCRF) on CodeProject. The HCRF is a discriminative model, forming the generative-discriminative pair with the hidden Markov model classifers.
This CodeProject article is a second on a series of articles about sequence classification, the first being about Hidden Markov Models. I've used this opportunity to write a little about generative versus discriminative models, and also provide a brief discussion on how Vapnik's ideas apply to these learning paradigms.
All the code available on those articles are also available within the Accord.NET Framework. Those articles provide good examples on how to use the framework and can be regarded as a practical implementation on how to use those models with the framework.
In the meantime, this article is also serving as a hook to a future article, an article about Hidden Conditional Random Fields (HCRFs). The HCRF models can serve the same purpose as the HMMs but can be generalized to arbitrary graph structures and be trained discriminatively, which could be an advantage on classification tasks.
As always, I hope readers can find it a good read :-)
With the beginning of this year, I would like to share a video I wish I had found earlier. It is about the recent breakthrough given by Deep Neural Networks in the field of speech recognition - which, despite I had known was a breakthrough, I didn't know it was already leading to such surprising great results.
Deep neural networks are also available in the Accord.NET Framework. However, they've been a very recent addition - if you find any issues, bugs, or just wish to collaborate on development, please let me know!
A new article has been published in CodeProject! This article details the Viola-Jones face detection algorithm available in the Accord.NET Framework. The article page also provides a standalone package for face detection which can be reused without instantiating the entire framework.
This month we have other great articles participating in the competition: Marcelo de Oliveira's Social News excels in the Web category; Roy's Inline MSIL in C# also presents a very interesting reading in the C# category. The later may eventually be extremely useful to leverage performance in managed applications.
For those who don't know, CodeProject is a amazingly useful site which publishes user-created articles and news. Every month, the best articles among all submissions are selected to win prizes and gain popularity as well!
I have manually translated and adapted the QuadProg solver for quadratic programming problems made by Berwin A. Turlach. His code was originally published under the GNU Library License, which has now been superseded by the GNU Lesser License. This adapted version honors the original work and is thus distributed under the same license.
Despite the name, the terms linear or quadratic programming have little resemblance to the set of activities most people now know as programming. Those terms usually usually refers to a specific set of function optimization methods, i.e. methods which can be used to determine the maximum or minimum points of special kinds of functions under a given number of solution constraints. For example, suppose we would like to determine the minimum value of the function:
f(x, y) = 2x + y + 4
Under the constraints that x and y must be non-negative (i.e. either positive or zero). This may seem fairly simple and trivial, but remember that practical linear programming problems may have hundreds or even thousands of variables and possibly million constraints.
When the problem to be solved involves a quadratic function instead of a linear function, but still presents linear constraints, this problem can be cast as a quadratic programming problem. Quadratic functions are polynomial functions in each each term may have at most a total degree of 2. For example, consider the function
f(x, y, z) = 2x² + 5xy + y² - z² + x – 5.
Now let's check the sum of the degrees for each variable on the polynomial terms. We start by writing the missing terms of the polynomial
and then proceed to check the sum of the degrees at each term. In the first term, 2+0+0 = 2. For the second, 1+1+0 = 2, and so on. Those functions have a nice property that they can be expressed in a matrix form
f(x) = 1/2 xT Qx + cTx.
Here, x and c are vectors. The matrix Q is a symmetric matrix specifying how the variables combine in the quadratic terms of the function. If this matrix is positive definite, then the function is convex, and the optimization has a single, unique optimum (we say it has a global optimum). The coefficients c specify the linear terms of our function.
Source code
The available source code is based on a translation of the Fortran code written by Berwin A. Turlach. However, some modifications have been made. Fortran uses column-major ordering for matrices, meaning that matrices are stored in memory in sequential order of column elements. Almost all other languages use row-major ordering, including C, C++, Java and C#. In order to improve data locality, I have modified the code to use the transpose of the original matrices D and A. I have also modified the QP formulation adopted in the Goldfarb and Idnani paper to reflect the form presented in the introduction.
This code is part of the Accord.NET Framework. However, the version available in this blog post will most likely not be the most recently, updated, fixed and enhanced version of the code. For the latest version, be sure to download the latest version of the framework on the project site or through a NuGet package.
Using the code
The first step in solving a quadratic programming problem is, well, specifying the problem. To specify a quadratic programming problem, one would need two components: a matrix D describing the relationship between the quadratic terms, and a vector d describing the linear terms. Perhaps this would work better with an example.
Suppose we are trying to solve a minimization problem. Given a function, the goal in such problems is to find the correct set of function arguments which would result in the minimum possible value for the function. An example of a quadratic minimization problem is given below:
However, note that this problem involves a set of constraints. The required solution for this minimization problem is required to lie in the interval specified by the constraints. More specifically, any x and y pair candidate for being a minimal of the function must respect the relations x - y = 5 and x >= 10. Thus, instead of lying in the unconstrained minimum of the function surface shown above, the solution lies slightly off the center of the surface. This is an obvious easy problem to solve manually, but it will fit for this demonstration.
As it can be seen (and also live demonstrated by asking Wolfram Alpha) the solution lies on the point (10,5), and the constrained minimum of the function is given by 170. So, now that we know what a quadratic programming problem looks like, how can we actually solve it?
Specifying the objective function
The first step in solving a quadratic programming problem is to specify the objective function. Using this code, there are three ways to specify it. Each of them has their own advantages and disadvantages.
1. Manually specifying the QP matrix.
This is the most common approach for numerical software, and probably the most cumbersome for the user. The problem matrix has to be specified manually. This matrix is sometimes denoted Q, D or H as it actually denotes the Hessian matrix for the problem.
The matrix Q is used to describe the quadratic terms of our problem. It is a n x n matrix, in which n corresponds to the number of variables in our problem, covering all possible combinations of variables. Recall our example given on the start of this section. We have 2 variables, x and y. Thus, our matrix Q is 2 x 2. The possible combinations for x and y are expressed in the table below.
x
y
x
x*x
x*y
y
y*x
y*y
To form our matrix Q, we can take all coefficients associated with each pair mentioned on the table above. The diagonal elements should also be multiplied by two (this is actually because the matrix is the Hessian matrix of the problem: it is the matrix of all second-order derivatives for the function. Since we have only at most quadratic terms, the elementary power rule of derivation “drops” the ² from the x² and y² terms – I think a mathematician would hit me with a stick for explaining it like this, but it serves well for a quick, non-technical explanation).
Remember our quadratic terms were 2x² - 1xy + 4y². Writing the terms on their proper position and differentiating, we have:
As it can be seen, the matrix is also symmetric (and often, but not always, positive definite). The next step, more trivial, is to write a vector d containing the linear terms. The linear terms are –5x –6y, and thus our vector d can be given by:
Therefore our C# code can be created like this:
double[,] Q = { { +4, -1 }, { -1, +8 }, };
double[] d = { -5, -6 };
2. Using lambda expressions
This approach is a bit more intuitive and less error prone. However, it involves lambdas functions and some people find it hard to follow them. Another disadvantage is that we will lose the edit & continue debugging ability of visual studio. The advantage is that the compiler may catch some obvious syntax errors automatically.
Note that the x and y variables could have been initialized to any value. They are only used as symbols, and not used in any computations.
3. Using text strings
This approach is more intuitive but a bit more error prone. The function can be specified using strings, as in a standard mathematical formula. However, since all we have are strings, there is no way to enforce static, compile time checking.
QuadraticObjectiveFunction f = new QuadraticObjectiveFunction("2x² - xy + 4y² - 5x - 6y");
Couldn’t be easier.
Specifying the constraints
The next step in specifying a quadratic programming problem is to specify the constraints. The constraints can be specified in almost the same way as the objective function.
1. Manually specifying the constraints matrix
The first option is to manually specify the constraints matrix A and vector b. The constraint matrix expresses the way the variables should be combined when compared to corresponding value on vector b. It is possible to specify either equality constraints or inequality constraints. The formulation used in this code is slightly different from the one used in Turlach’s original implementation. The constraints are supposed to be in the form:
A1x = b1
A2x = b2
This means that each line of matrix A expresses a possible combination of variables x which should be compared to the corresponding line of b. An integer variable m can be specified to denote how many of the first rows of matrix A should be treated as equalities rather than inequalities. Recall that in our example the constraints are given by 1x -1y = 5 and 1x = 10. Lets write this down in a tabular form:
#
x
y
?
b
q1
1
-1
=
5
q2
1
0
=
10
Thus our matrix A and vector b can be specified as:
And not forgetting that m = 1, because the first constraint is actually an equality.
2. Using classes and objects
A more natural way to specify constraints is using the classes and objects of the Accord.NET Framework. The LinearConstraint class allows one to specify a single constraint using an object-oriented approach. It doesn’t have the most intuitive usage on earth, but has much more expressiveness. It can also be read aloud, it that adds anything! :-)
List<LinearConstraint> list = new List<LinearConstraint>();
list.Add(new LinearConstraint(numberOfVariables: 2) { VariablesAtIndices = newint[] { 0, 1 }, // index 0 (x) and index 1 (y) CombinedAs = newdouble[] { 1, -1 }, // when combined as 1x -1y ShouldBe = ConstraintType.EqualTo, Value = 5 });
The specification is centered around the notion that variables are numbered and have an associated index. For example, x is the zero-th variable of the problem. Thus x has an index of 0 and y has an index of 1. So for example, reading aloud the last constraint, it is possible to express how the variables at indices 0 and 1, when combined as 1x and –1y, should be equal to value 5.
2. Using lambda expressions
A more intuitive way to express constraints is again using lambda expressions. And again the problems are the same: some people find it hard to follow and we lose edit & continue.
var constraints = new List<LinearConstraint>(); constraints.Add(new LinearConstraint(f, () => x - y == 5)); constraints.Add(new LinearConstraint(f, () => x >= 10));
3. Using text strings
Same as above, but with strings.
var constraints = new List<LinearConstraint>(); constraints.Add(new LinearConstraint(f, "x - y = 5")); constraints.Add(new LinearConstraint(f, "x >= 10"));
Finally, creating and solving the problem
Once we have specified what do we want, we can now ask the code for a solution. In case we have opted for manually specifying the matrix A, vector b and integer m, we can use:
// Create the optimization problem var solver = new GoldfarbIdnaniQuadraticSolver(numberOfVariables:2, A, b, m);
In case we have opted for creating a list of constraints instead, we can use:
// Create our optimization problem var solver = new GoldfarbIdnaniQuadraticSolver(numberOfVariables: 2, constraints: list);
After the solver object has been created, we can call Minimize() to solve the problem. In case we have opted for manually specifying Q and d, we can use:
// Attempt to solve the problem double minimumValue = solver.Minimize(Q, d);
And in case we have opted for creating a QuadraticObjectiveFunction object, we can use:
// Attempt to solve the problem double minimumValue = target.Minimize(f);
In either case, the solution will be available in the Solution property of the solver object, and will be given by:
double value = solver.Value; // f(x,y) = 170 double x = solver.Solution[0]; // x = 10 double y = solver.Solution[1]; // y = 5
Because the code has been translated by hand (in contrast of using automatic translators such as f2c) there could be potential bugs in the code. I have tested the code behavior against R’s quadprog package and still didn’t find errors. But this does not mean the code is bug-free. As always, as is the case of everything else in this blog, this code is published in good faith, but I can not guarantee the correctness of everything. Please read the disclaimer for more information.
D. Goldfarb and A. Idnani. Dual and Primal-Dual Methods for Solving Strictly Convex Quadratic Programs (1982). In J. P. Hennart (ed.), Numerical Analysis, Springer-Verlag, Berlin, pages 226–239.
A new version (2.2.0) of the Accord.NET Framework has just been released. This new version introduces many new features, fixes and improvements. The most interesting additions are certainly the HeadController and FaceController .NET components.
Accord.NET Framework sample application for Gesture Controller Components
The Accord.NET Controller components can be used to generate events based on webcam motion. By using a combination of HaarCascadeClassifiers, Camshift and Template-based Tracking, those components are able to detect when a face enters scene, leaves the scene, and moves across a scene.
The video above shows only the sample application which comes together with the framework. However, the interesting part is that this is just a sample of what can be accomplished using the real controller components. The controller components are .NET components, similar to Button, Label or Timer, and can be dragged and dropped from Visual Studio’s ToolBox directly into any application.
Once inside an application, it will be possible to set event actions just as in any other .NET component:
The controls have built-in support for calibration. All values except tilting angle are passed to the hosting application in the [-1;+1] range, in which -1 indicates either a total left/down/backwards position and +1 indicates a total right/up/forward position. The tilting angle is given in radians. Please note that the face controller is still a bit experimental and still requires some tuning.
This new version also introduces HSL Color Range object trackers, more default Haar Cascades, an experimental version of linear-chain Conditional Random Fields, and the ability to generate hardcoded C# definitions of any Haar cascade available in the OpenCV XML format. There is also initial support for finger detection using new implementations for Border-Following contour extraction, K-Curvatures and Convex Hull Defects extraction. On the statistics side, there has been the inclusion of the Von-Mises distribution, Moving and Running Normal distributions and improvements in the Multivariate Gaussian implementation. The full release notes are available in the release's download page.
Also, a special thanks to Professor Dr. Modesto Castrillón Santana to let me embed some of his Haar definitions into the framework under the LGPL license. Please be sure to include a reference to his work if you plan to use this in an academic publication.
As always, I hope those additions and improvements will be useful to everyone :-)
The next version of the Accord.NET Framework will feature two important additions, alongside with a new namespace to accommodate them: The Camshift object tracker and the Viola-Jones object detector. Both will be located inside the new Accord.Vision namespace for Computer Vision algorithms.
Camshift object tracker
Viola-Jones object detector
Additionally, other cool additions include the availability of Multi-class Kernel Support Vector Machines (with support for parallelized learning algorithms), Generic Cross-validation classes for model evaluation and Generic Gridsearch classes for model parameter tuning.
The Linear and Non-linear Kernel Discriminant Analysis implementations have been updated to use the Generalized Eigenvalue Decomposition instead of the Eigenvalue Decomposition following a prior matrix inversion. This has cut execution time in nearly half. Other optimizations have also been made to some radial basis Kernel functions, although they are not very noticeable.
The next version of the Accord.NET Framework will be labeled version 2.1.0. There are some interface changes that may break compatibility with older versions. Also the minimum required version of AForge.NET Framework has been leveled up to 2.1.3. You may have to update your assemblies if you wish to upgrade an already existing project.
Today I have just released the first version of Accord.NET. The Accord.NET Framework is a C# framework I have developed over the years while working on several areas of artificial intelligence, machine learning and statistics.
The Accord.NET Framework extends the excellent AForge.NET Framework with new tools and libraries. In order to use Accord.NET, you must have AForge.NET already installed. The first version number of Accord.NET will be used to indicate the compatibility status with AForge.NET versions, thus the first version will be starting at 2.0.0.
The framework is comprised by the set of libraries and sample applications, which demonstrate their features:
Accord.Statistics - library with statistical analysis and other tools;
Accord.Imaging - extension to the AForge.NET Imaging library with new filters and routines;
Accord.Neuro - extension to the AForge.NET Neuro library with other learning algorithms;
Accord.MachineLearning - extension to AForge's machine learning library with Support Vector Machines;
Accord.Audio - experimental library with filters and audio processing routines.
The Framework has just been released, so be ready to expect bugs and unpolished/unfinished areas. The code is released under a LGPL license. For additional help, issues, and discussions, please refer to the recently created forums.
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.
The source code available here is distributed under a GPL license. The sequential minimal optimization for regression was implemented by following Smola’s Tutorial on Support Vector Regression. However, modifications had been based on GPL code by Sylvain Roy.
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 z1, z2, … , zN and weights w1, w2, … , wN is then given by:
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 z1, z2, … , zN and weights w1, w2, … , wN is now given by:
It is also very straightforward to see that, using a linear kernel of the form K(z,x) = <z,x> = zTx, we recover the original formulation for the linear SVM.
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.
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.
Polynomial Kernel
The Polynomial kernel is a non-stationary kernel. It is well suited for problems where all data is normalized.
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.
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.
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'sSequential 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 theother 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 yi can be individually multiplied by their corresponding Lagrange multipliers ai to form a single weight vector w.
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 ai 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 ai are then subtracted to form a single weight vector w.
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.
/// <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> publicint Inputs { get { return inputCount; } }
/// <summary> /// Gets or sets the collection of support vectors used by this machine. /// </summary> publicdouble[][] SupportVectors { get { return supportVectors; } set { supportVectors = value; } }
/// <summary> /// Gets or sets the collection of weights used by this machine. /// </summary> publicdouble[] Weights { get { return weights; } set { weights = value; } }
/// <summary> /// Gets or sets the threshold (bias) term for this machine. /// </summary> publicdouble 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> publicvirtualdouble 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> publicdouble[] Compute(double[][] inputs) { double[] outputs = newdouble[inputs.Length];
for (int i = 0; i < inputs.Length; i++) outputs[i] = Compute(inputs[i]);
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] publicclass 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) thrownew 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> /// publicoverridedouble Compute(double[] inputs) { double s = Threshold;
for (int i = 0; i < SupportVectors.Length; i++) s += Weights[i] * kernel.Function(SupportVectors[i], inputs);
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> /// publicclass SequentialMinimalOptimization : ISupportVectorMachineLearning { privatestatic Random random = new Random();
// Training data privatedouble[][] inputs; privateint[] outputs;
// Error cache to speed up computations privatedouble[] 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) thrownew ArgumentNullException("machine");
if (inputs == null) thrownew ArgumentNullException("inputs");
if (outputs == null) thrownew ArgumentNullException("outputs");
if (inputs.Length != outputs.Length) thrownew 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) thrownew 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) thrownew 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> publicdouble Complexity { get { returnthis.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> publicbool 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> publicdouble 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> publicdouble Tolerance { get { returnthis.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> /// publicdouble 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();
// 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); }
// 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 = newdouble[vectors][]; machine.Weights = newdouble[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;
/// <summary> /// Runs the SMO algorithm. /// </summary> /// /// <returns> /// The misclassification error rate of /// the resulting support vector machine. /// </returns> /// publicdouble Run() { return Run(true); }
/// <summary> /// Computes the error rate for a given set of input and outputs. /// </summary> /// publicdouble 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> /// privateint 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> /// privatebool takeStep(int i1, int i2) { if (i1 == i2) returnfalse;
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); }
/// <summary> /// Computes the SVM output for a given point. /// </summary> /// privatedouble 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; }
privatedouble 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; }
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.
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);
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 SequentialMinimalOptimization learn = new SequentialMinimalOptimization( machine, inputs, xor);
// Run the learning algorithm learn.Run();
After the model has been trained, we can compute its outputs for the given inputs.
double[] output = machine.Compute(inputs);
The machine should be able to correctly identify all of the input instances.
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.
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.
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.
.png)
