For quite some time I have been trying to find a Java or C# implementation of the Generalized Eigenvalue Decomposition. I just wanted something simple, that didn’t depend on any external libraries such as BLAS or LAPACK. As I could not find any simple and understandable implementation of this decomposition, I decided to implement one myself.
where A and B are matrices. If v obeys this equation, with some λ, then we call v the generalized eigenvector of A and B, and λ is called the generalized eigenvalue of A and B which corresponds to the generalized eigenvector v. The possible values of λ must obey the following equation
An interesting feature of the generalized eigenvalue decomposition is that it finds the eigenvectors of the matrix B-1A even if B is singular and does not have an inverse. If B is nonsingular, the problem could be solved by reducing it to a standard eigenvalue problem of the form B-1Ax=λx. However, because B can be singular, an alternative algorithm, called the QZ method, is necessary.
EISPACK is a software library for numerical computation of eigenvalues and eigenvectors of matrices, written in FORTRAN. It was originally written around 1972–1973, based heavily on algorithms originally implemented in ALGOL. Because those algorithms are much more human-readable than their LAPACK counterparts, I decided to port the related functions from EISPACK.
The functions which implement the Generalized Eigenvalue Decomposition in EISPACK are called QZHES, QZIT, QZVAL and QZVEC. As their name implies, they use the QZ method for finding the generalized eigenvalues of a matrix pair (A,B).
Source code
The source code is available in the download link in the upper part of this article. It can also be found in the latest versions of the Accord.NET Framework. The algorithm presented in the code is equivalent to the eig(A,B) command of MATLAB.
Because the code has been ported from FORTRAN, many sections of the algorithms have stayed as they appear in their original functions. Please be aware that it may contain lots of labels and goto's, as most structured programming features we take for granted today were not commonplace at the time.
Using the code
Using the code is much simpler than dealing with LAPACK functions. Because the functions have been implemented directly in C#, the code is somewhat more human readable and has also been wrapped in the same nice object-oriented approach as the other matrix decompositions that can be found in MAPACK and JAMA.
Using the given sources, the following MATLAB code:
[V,D] = eig(A,B)
Becomes equivalent to the C# code:
var gevd = new GeneralizedEigenvalueDecomposition(A,B); var V = gevd.Eigenvectors; var D = gevd.DiagonalMatrix;
in the sense that both satisfy the identity A*V = B*V*D .
RANSAC is an iterative method to build robust estimates for parameters of a mathematical model from a set of observed data which is known to contain outliers. The RANSAC algorithm is often used in computer vision, e.g., to simultaneously solve the correspondence problem and estimate the fundamental matrix related to a pair of stereo cameras.
RANSAC is an abbreviation for "RANdom SAmple Consensus". It is an iterative method to estimate parameters of a mathematical model from a set of observed data which may contains outliers. It is a non-deterministic algorithm in the sense that it produces a reasonable result only with a certain probability, with this probability increasing as more iterations are allowed. The algorithm was first published by Fischler and Bolles in 1981.
The basic assumption is that the data consists of "inliers", i.e., data whose distribution can be explained by some mathematical model, and "outliers" which are data that do not fit the model. Outliers could be considered points which come from noise, erroneous measurements or simply incorrect data. RANSAC also assumes that, given a set of inliers, there exists a procedure which can estimate the parameters of a model that optimally explains or fits this data.
Example: Fitting a simple linear regression
We can use RANSAC to robustly fit a linear regression model using noisy data. Consider the example below, in which we have a cloud of points that seems to belong to a line. These are the inliers of the data. The other points, which can be seem as measurement errors or extreme noise values, are points expected to be considered outliers.
Linear structure contained in noisy data.
RANSAC is able to automatically distinguish the inliers from the outliers through the evaluation of the linear regression model. To do so, it randomly selects subsets from the data and attempts to fit linear regression models using them. The model which best explains most of the data will then be returned as the most probably correct model fit.
The image below shows the result of fitting a linear regression directly (as shown by the red line) and using RANSAC (as shown by the blue line). We can see that the red line represents poorly the data structure because it considers all points in order to fit the regression model. The blue line seems to be a much better representation of the linear relationship hidden inside the overall noisy data.
Hidden linear structure inside noisy data. The red line shows the fitting of a linear regression model directly considering all data points. The blue line shows the same result using RANSAC.
Source code
The code below implements RANSAC using a generic approach. Models are considered to be of the reference type TModel and the type of data manipulated by this model is considered to be of the type TData. This approach allows for the creation of a general purpose RANSAC algorithm which can be used in very different contexts, be it the fitting of linear regression models or the estimation of homography matrices from pair of points in different images.
/// <summary> /// Computes the model using the RANSAC algorithm. /// </summary> public TModel Compute(TData[] points, outint[] inliers) { // We are going to find the best model (which fits // the maximum number of inlier points as possible). TModel bestModel = null; int[] bestInliers = null; int maxInliers = 0;
// For this we are going to search for random samples // of the original points which contains no outliers.
int count = 0; // Total number of trials performed double N = maxEvaluations; // Estimative of number of trials needed.
// While the number of trials is less than our estimative, // and we have not surpassed the maximum number of trials while (count < N && count < maxEvaluations) { int[] idx; TModel model = null; int samplings = 0;
// While the number of samples attempted is less // than the maximum limit of attempts while (samplings < maxSamplings) { // Select at random s datapoints to form a trial model. idx = Statistics.Tools.Random(points.Length, s); TData[] sample = points.Submatrix(idx);
// If the sampled points are not in a degenerate configuration, if (!degenerate(sample)) { // Fit model using the random selection of points model = fitting(sample); break; // Exit the while loop. }
samplings++; // Increase the samplings counter }
// Now, evaluate the distances between total points and the model returning the // indices of the points that are inliers (according to a distance threshold t). idx = distances(model, points, t);
// Check if the model was the model which highest number of inliers: if (idx.Length > maxInliers) { // Yes, this model has the highest number of inliers.
maxInliers = idx.Length; // Set the new maximum, bestModel = model; // This is the best model found so far, bestInliers = idx; // Store the indices of the current inliers.
// Update estimate of N, the number of trials to ensure we pick, // with probability p, a data set with no outliers. double pInlier = (double)idx.Length / (double)points.Length; double pNoOutliers = 1.0 - System.Math.Pow(pInlier, s);
N = System.Math.Log(1.0 - probability) / System.Math.Log(pNoOutliers); }
count++; // Increase the trial counter. }
inliers = bestInliers; return bestModel; }
Besides the generic parameters, the class utilizes three delegated functions during execution.
The Fitting function, which should accept a subset of the data and use it to fit a model of the chosen type, which should be returned by the function;
The Degenerate function, which should check if a subset of the training data is already known to result in a poor model, to avoid unnecessary computations; and
The Distance function, which should accept a model and a subset of the training data to compute the distance between the model prediction and the expected value for a given point. It should return the indices of the points only whose predicted and expected values are within a given threshold of tolerance apart.
Using the code
In the following example, we will fit a simple linear regression of the form x→y using RANSAC. The first step is to create a RANSAC algorithm passing the generic type parameters of the model to be build, i.e. SimpleLinearRegression and of the data to be fitted, i.e. a double array.
In this case we will be using a double array because the first position will hold the values for the input variable x. The second position will be holding the values for the output variables y. If you are already using .NET 4 it is possible to use the Tuple type instead.
// Create a RANSAC algorithm to fit a simple linear regression var ransac = new RANSAC<SimpleLinearRegression, double[]>(minSamples); ransac.Probability = probability; ransac.Threshold = errorThreshold; ransac.MaxEvaluations = maxTrials;
After the creation of the RANSAC algorithm, we should set the delegate functions which will tell RANSAC how to fit a model, how to tell if a set of samples is degenerate and how to check for inliers in data.
// Set the RANSAC functions to evaluate and test the model
ransac.Fitting = // Define a fitting function delegate(double[][] sample) { // Retrieve the training data double[] inputs = sample.GetColumn(0); double[] outputs = sample.GetColumn(1);
// Build a Simple Linear Regression model var r = new SimpleLinearRegression(); r.Regress(inputs, outputs); return r; };
ransac.Degenerate = // Define a check for degenerate samples delegate(double[][] sample) { // In this case, we will not be performing such checkings. returnfalse; };
ransac.Distances = // Define a inlier detector function delegate(SimpleLinearRegression r, double[][] sample, double threshold) { List<int> inliers = new List<int>(); for (int i = 0; i < sample.Length; i++) { // Compute error for each point double input = sample[i][0]; double output = sample[i][1]; double error = r.Compute(input) - output;
// If the squared error is below the given threshold, // the point is considered to be an inlier. if (error * error < threshold) inliers.Add(i); } return inliers.ToArray(); };
Finally, all we have to do is call the Compute method passing the data. The best model found will be returned by the function, while the given set of inliers indices for this model will be returned as an out parameter.
// Finally, try to fit the regression model using RANSAC int[] idx; SimpleLinearRegression rlr = ransac.Compute(data, out idx);
Sample application
The accompanying source application demonstrates the fitting of the simple linear regression model with and without using RANSAC. The application accepts Excel worksheets containing the independent values in the first column and the dependent variables in the second column.
References
Wikipedia contributors. RANSAC. Wikipedia, The Free Encyclopedia.
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.
One of the firsts operators for interest point detection was developed by Hans P. Moravec in 1977 for his research involving the automatic navigation of a robot through a clustered environment. It was also Moravec who defined the concept of "points of interest" in a image and concluded these interest points could be used to find matching regions in different images.
The Moravec operator is considered to be a corner detector because it defines interest points as points where there are large intensity variations in all directions. This often is the case at corners. It is interesting to note, however, that Moravec was not specifically interested in finding corners, just distinct regions in an image that could be used to register consecutive image frames.
Harris and Stephens improved upon Moravec's corner detector by considering the differential of the corner score with respect to direction directly. The Harris corner detector computes the locally averaged moment matrix computed from the image gradients, and then combines the Eigenvalues of the moment matrix to compute a corner measure, from which maximum values indicate corners positions.
Below is the source code for the Harris Corners Detector algorithm. The code utilizes the excellent AForge.NET Framework and is also molded in its same standards.
/// <summary> /// Process image looking for corners. /// </summary> /// /// <param name="image">Source image data to process.</param> /// /// <returns>Returns list of found corners (X-Y coordinates).</returns> /// /// <exception cref="UnsupportedImageFormatException"> /// The source image has incorrect pixel format. /// </exception> /// publicunsafe List<IntPoint> ProcessImage(UnmanagedImage image) { // check image format if ( (image.PixelFormat != PixelFormat.Format8bppIndexed) && (image.PixelFormat != PixelFormat.Format24bppRgb) && (image.PixelFormat != PixelFormat.Format32bppRgb) && (image.PixelFormat != PixelFormat.Format32bppArgb) ) { thrownew UnsupportedImageFormatException("Unsupported pixel format of the source image."); }
// make sure we have grayscale image UnmanagedImage grayImage = null;
// get source image size int width = grayImage.Width; int height = grayImage.Height; int srcStride = grayImage.Stride; int srcOffset = srcStride - width;
fixed (float* pdx = diffx, pdy = diffy, pdxy = diffxy, pmap = map) { float* dx = pdx; float* dy = pdy; float* dxy = pdxy; float* H = pmap; float M, A, B, C;
for (int y = 0; y < height; y++) { for (int x = 0; x < width; x++, dx++, dy++, dxy++, H++) { A = *dx; B = *dy; C = *dxy;
if (measure == HarrisCornerMeasure.Harris) { // Original Harris corner measure M = (A * B - C * C) - (k * ((A + B) * (A + B))); } else { // Harris-Noble corner measure M = (A * B - C * C) / (A + B + Accord.Math.Special.SingleEpsilon); }
if (M > threshold) { *H = M; // insert value in the map } } } }
// 4. Suppress non-maximum points List<IntPoint> cornersList = new List<IntPoint>();
// for each row for (int y = r, maxY = height - r; y < maxY; y++) { // for each pixel for (int x = r, maxX = width - r; x < maxX; x++) { float currentValue = map[y, x];
// for each windows' row for (int i = -r; (currentValue != 0) && (i <= r); i++) { // for each windows' pixel for (int j = -r; j <= r; j++) { if (map[y + i, x + j] > currentValue) { currentValue = 0; break; } } }
// check if this point is really interesting if (currentValue != 0) { cornersList.Add(new IntPoint(x, y)); } } }
return cornersList; }
Using the code
Code usage is very simple. A Corners Marker filter from the AForge Framework can be used to directly draw the detected interest points in the original image as it would be usual within the framework.
// Open a image Bitmap image = ...
// Create a new Harris Corners Detector using the given parameters HarrisCornersDetector harris = new HarrisCornersDetector(k); harris.Threshold = threshold; harris.Sigma = sigma;
// Create a new AForge's Corner Marker Filter CornersMarker corners = new CornersMarker(harris, Color.White);
// Apply the filter and display it on a picturebox pictureBox1.Image = corners.Apply(image);
Sample application
The accompanying sample application is pretty much self-explanatory. It performs the corners detection in the famous Lena Söderberg's picture, but can be adapted to work with other pictures as well. Just add another image to the project settings' resources and change the line of code which loads the bitmap.
It is also very straightforward to adapt the application to load arbitrary images from the hard disk. I have opted to leave it this way for simplicity.
The Code Project is a free repository of source code, tutorials and development resources. It is also home of a large and ever-growing community of software developers, architects and other technology professionals. It has fairly active forums, and is a reasonably good resource for resolving difficult software development issues.
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.
