Matrix manipulation using Accord.NET

Matrix manipulation in Accord.NET
is very straightforward: Just add a new using directive on top of your class to
have (literally) about a hundred new extension methods that operate directly on
.NET multi-dimensional arrays.

target="_blank">
width="358" height="265" />

Introduction

target="_blank">
border="0" alt="accord-matrix" align="right" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg5TMxjbp7bYPTPK7yI2QZyk-40mRVxet3XhWaSWb2Udc9n6ASMrz13p7rCDrQTeXpHdY_jZnehvXaPjjYbSdhlAI29tgvIYLqTc9iTHrH4sKAUoE-qUIamPv_itkZgo2Y1mx1UTXW6mkhZ/?imgmax=800"
width="246" height="251" />

Accord.NET uses a bit different approach for matrix manipulation in contrast to
other libraries. By using target="_blank">C# 3.0 extension methods, Accord adds several of the standard
methods you would expect from a Matrix library,
such as linear system solving, matrix algebra and numerical decompositions directly
to the standard double[,] (or more generally T[,]) matrices of the framework[ href="http://accord-net.origo.ethz.ch/image/class_diagram_accord_math">^].

This approach offers some advantages since it avoids mismatches when one is using
multiple libraries, each with their own and often incompatible Matrix implementations.
Extending multi-dimensional arrays makes the use of matrices much more natural in
the .NET world, as no specialized Matrix classes have to be used and no code has
to be rewritten just to use a particular Matrix implementation.

 

Please note, however, that most methods implemented by Accord are not equivalent
to the heavily optimized versions of more specialized numerical packages, such as

BLAS, from a performance view. Their real power comes when prototyping or realizing algorithms into code. Once an algorithm is written, several
unit tests can be created to test the correctness of the code. After that, because
we have an early working prototype, it will be much easier to perform optimizations
as needed. The key point is to have a working version early which can (if
actually necessary) be optimized later.

Using extension methods

The first step is to include a new using directive on the top of
your source file.

title="accord-matrix3" border="0" alt="accord-matrix3" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgpo_b_J445BwyGME6LhJeH5ub3pRV4uJ5zYb0sSgXhXXdqC4m6HkzlvYv_8ZCEfJc_RssdRmbM_NbgE245V4a8yZgMu9Xohx8R6cEdKbO4D-7GRVciWE3fzE7_Kx4TBzg7llTX1PaZkzT0/?imgmax=800"
width="377" height="137" />

This is all that is necessary after you have referenced the Accord.Math library
in your project. Now you can use all of the methods and operations described below.
The following list is nowhere complete, but shows only some basics and highlights
from the current version of Accord.

Declaring matrices

Using standard .NET declarations

To declare matrices, no special code is required. Just create multi-dimensional
arrays as usual.

double[,] A = { {1, 2, 3}, {6, 2, 0}, {0, 0, 1} };

double[,] B = { {2, 0, 0}, {0, 2, 0}, {0, 0, 2} };

Using Accord.NET extensions

Additionally, one may wish to use one of the convenience methods of Accord.NET to
create specialized matrices such as the Identity matrix, multiples of the Identity
matrix or Diagonal matrices.

Accord.NET (C#)	MATLAB®
double[,] B = class="str">Matrix.Identity(3).Multiply(2);	B = eye(3)*2;
double[,] B = class="str">Matrix.Diagonal(3, 2.0);	B = eye(3)*2;

Using Accord.NET extensions with implicit typing

By using implicit type variable declaration, the code acquires a certain lightweight
feeling and gets closer of its MATLAB®/Octave counterpart (MATLAB is a registered
trademark of The MathWorks, Inc).

Accord.NET (C#)	MATLAB®
var I = class="str">Matrix.Identity(3);	I = eye(3)

A much more closer alternative will be possible by using Algorithm Environments,
a upcoming feature in Accord. When available, this feature will allow construction
of mathematical code using

var I =         class="kwrd">eye(3)

directly. The code will be based on current Octave syntax. Octave is a high-level
language, primarily intended for
numerical computations, that is mostly compatible with
MATLAB.

Matrix operations

All standard matrix operations such as transpose, inverse, column and row manipulations
are available in the extension methods. In the example below, A is the same standard
double[,] matrix declared in the first section of this
article. All methods return a new double[,] matrix
as a result, leaving the original matrix A untouched.

Operation	Accord.NET (C#)	MATLAB®
Transpose	var At = A.Transpose();	At = A'
Inverse	var invA = A.Inverse();	invA = inv(A)
Pseudo-Inverse	var pinvA = A.PseudoInverse();	pinvA = pinv(A)

Matrix Algebra

All common algebraic operations for matrices are also implemented. Those are the
common operations such as addition, subtraction and multiplication. Here, division
is actually a shortcut for multiplying by the inverse.

Operation	Accord.NET (C#)	MATLAB®
Addition	var C = A.Add(B);	C = A + B
Subtraction	var C = A.Subtract(B);	C = A - B
Multiplication	var C = A.Multiply(B);	C = A * B
Division	var C = A.Divide(B);	C = A / B

The above also works with vectors and scalars.

Operations	Accord.NET (C#)	MATLAB®
Multiplying by a scalar	var H = A.Multiply(3.14);	H = A * 3.14;
Multiplying by a column vector	double[] u = { 1, 2, 3 }; double[] v = A.Multiply(u);	v = A * u';

Special element-wise operations

Element-wise operations are operations performed on each element of the matrix or
vector. In Matlab, they are some times known as the dot operations, since they are
usually denoted by prefixing a dot on the common operators.

Operation	Accord.NET (C#)	MATLAB®
Multiplication	var C = A.ElementwiseMultiply(B);	C = A .* B
Division	var C = A.ElementwiseDivide(B);	C = A ./ B
Power	var C = A.ElementwisePower(B);	C = A .^ B

Vector operations

Accord.NET can also perform many vector operations. Among them are the many flavors
of products between vectors, such as the target="_blank">inner, the target="_blank">outer and the target="_blank">Cartesian.

Operation	Accord.NET (C#)	MATLAB®
Inner product (a.k.a. the scalar product)	var w = u.InnerProduct(v);	w = u*v'
Outer product (a.k.a. the matrix product)	var w = u.OuterProduct(v);	w = u'*v
Cartesian product	var w = u.CartesianProduct(v);
Euclidean Norm	double n = u.Norm();	n = norm(u)
Sum	double s = u.Sum();	s = sum(u)
Product	double p = u.Product();	p = prod(u)

Matrix characteristics

Some common matrix characteristics, such as the determinant and trace, are readily
available.

Operation	Accord.NET (C#)	MATLAB®
Determinant	A.Determinant();	det(A)
Trace	A.Trace();	tr(A)

Other characteristics

Other available characteristics are the Summation and Product of vector and matrix
elements.

Operation	Accord.NET (C#)	MATLAB®
Sum vector	double[] sum = A.Sum();	sum(A)
Sum of elements	double sum = A.Sum().Sum()	sum( class="kwrd">sum(A))
Sum along columns	double[] sum = A.Sum(0)	sum(A, 1)
Sum along rows	double[] sum = A.Sum(1)	sum(A, 2)

Linear Algebra

Linear algebra
is certainly one of the most important fields of mathematics, if not the most important
one. Accord includes many methods for numerical linear algebra such as matrix inversion
and matrix decompositions. Most of them were originally based on target="_blank">JAMA and MAPACK, but today Accord has some additions from
routines translated from EISPACK (mainly the
Generalized Eigenvalue Decomposition[^],
which is currently absent from Jama).

Operation	Accord.NET (C#)	MATLAB®
Solve a linear system	x = A.Solve(B)	x = A \ B

Eigenvalue Decomposition

Operation	Accord.NET (C#)	MATLAB®
Standard decomposition	var evd = new EigenvalueDecomposition(A); var V = evd.Eigenvectors; var D = evd.DiagonalMatrix;	[V, D] = eig(A)
Generalized decomposition	var gevd = new GeneralizedEigenvalueDecomposition(A,B); var V = gevd.Eigenvectors; var D = gevd.DiagonalMatrix;	[V, D] = eig(A,B)

Singular Value Decomposition

Operation	Accord.NET (C#)	MATLAB®
Economy decomposition	var svd = new SingularValueDecomposition(A); var U = evd.LeftSingularVectors; var S = evd.Diagonal; var V = evd.RightSingularVectors;	[U,S,V] = svd(A,0)

QR Decomposition

Operation	Accord.NET (C#)	MATLAB®
Standard decomposition	var qr = new QrDecomposition(A); var Q = qr.OrthogonalFactor; var R = qr.UpperTriangularFactor;	[Q,R] = QR(A)

Cholesky decomposition

Operation	Accord.NET (C#)	MATLAB®
Standard decomposition	var chol = new CholeskyDecomposition(A); var R = chol.LeftTriangularFactor;	R = CHOL(A)

LU Decomposition

Operation	Accord.NET (C#)	MATLAB®
Standard decomposition	var lu = new LuDecomposition(A); var L = lu.LowerTriangularFactor; var U = lu.UpperTriangularFactor;	[L, U] = LU(A)

Special operators

There are some other useful operators available in Accord.NET. There are facilities
to create index vectors (very common in Matlab for accessing sub-portions of a matrix),
select elements, find elements based on a selection criteria, and so on.

Operation	Accord.NET (C#)	MATLAB®
Create a vector of indices	var idx = class="str">Matrix.Indices(0,10);	idx = 1:9
Selecting elements	var B = A.Submatrix(idx);	B = A(idx)
Finding elements matching a certain criteria (For example, finding x ∈ v / x > 2).	var v = { 5, 2, 2, 7, 1, 0 }; var idx = v.Find(x => x > 2);	v = [ 5 2 2 7 1 0]; idx = find(v > 2)
Reshaping a vector to a matrix.	double[] m = { 1, 2, 3, 4 }; double[,] M = Matrix.Reshape(m, 2, 2);	m = [1 2 3 4]; reshape(m,2,2)

 

More than just matrices

Accord.NET also offers many other standard mathematical functions such as the Gamma
function, log-Gamma, Digamma, Bessel functions, Incomplete beta integrals, and so
on.

target="_blank">
title="accord-matrix4" border="0" alt="accord-matrix4" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgBD_5qI0IflDtm_E5eejzkOOLYsMwGu5OhTlsJ4YNqLJvQFmYNtIG0LocelOaRoNKe0iLAyny0EaKwiYQ_VEVMg-KGXXCWllFpM3XXpQDDHFuF739_mvxNbb1e2sF34hId6sEIraixKkOh/?imgmax=800"
width="412" height="357" />

For a complete listing of the framework features, please check the target="_blank">project page at Origo. In particular, don’t forget to check
out the the target="_blank">class diagram for the Accord.Math namespace.

Cheers!

 

Legal notice: MATLAB is a registered trademark of The MathWorks, Inc. All
other trademarks are the property of their respective owners.

Automatic Image Stitching with Accord.NET

I have just posted a new article on CodeProject, entitled Automatic Image Stitching with Accord.NET. It is a demonstration of automatic image stitching by interest point matching using the Accord and AForge.NET Frameworks.

 

 

The method used is pretty classic in the computer vision literature. There are many other variations of the method with their own advantages and disadvantages, but most of them are built around the same ideas - feature detection, matching and blending. This is one of the most straightforward and free implementations available, since some of those most sophisticated methods are patented.

New additions to Accord.NET: Computer Vision namespace, Camshift and Viola-Jones Detector

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.

Generalized Eigenvalue Decomposition in C#

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.

• Download source code.

Introduction

As Wikipedia says, a generalized Eigen value problem is the problem of finding a vector v that obeys

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.

Solving the problem using the QZ method is equivalent to computing the Eigen decomposition of B^-1A without the need of inverting B, which could be impossible or ill-conditioned if B is singular or near-singular. It is also computationally simpler as it saves the time and memory needed for inverting and storing the inverse of B.

EISPACK

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 .

References

• 
  

  Wolfram Research, “Generalized Eigenvalue problem”, Mathematica Documentation Center, 2010.

  
  


• 
  

  Wikipedia contributors, "Eigendecomposition of a matrix," Wikipedia, The Free Encyclopedia, 2010.

  
  

RANdom Sample Consensus (RANSAC) in C#

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.

• Download source code


• Download sample application


• Browse online the application source code or see also the framework's inner workings.

This code has also been incorporated in Accord.NET Framework, which includes the latest version of this code plus many other statistics and machine learning tools.

Introduction

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, out int[] 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.
            return false;
        };

    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.




• P. D. Kovesi. MATLAB and Octave Functions for Computer Vision and Image Processing. School of Computer Science & Software Engineering, The University of Western Australia.

Accord.NET Framework - An extension to AForge.NET

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.

Science!