Mathematics

Matrices

Declaring and using matrices in the Accord.NET Framework does not requires much. In fact, it does not require anything else that is not already present at the .NET Framework. If you have already existing and working code using other libraries, you don't have to convert your matrices to any special format used by Accord.NET. This is because Accord.NET is built to interoperate with other libraries and existing solutions, relying solely on default .NET structures to work.

To begin, please add the following using directive on top of your .cs (or equivalent) source code file:

using Accord.Math;

After this, all common .NET datatypes will be extended with several extension methods related to mathematics.

http://accord.googlecode.com/svn/wiki/samples/accord-math-matrixoperations-img.png

Creating matrices

The most straightforward way to declare a matrix in Accord.NET is simply using:

 double[,] matrix = 
 {
    { 1, 2 },
    { 3, 4 },
    { 5, 6 },
};

Yep, that is right. You don't need to create any fancy custom Matrix classes or vectors to make Accord.NET work, which is a plus if you have already existent code using other libraries. You are also free to use both the multidimensional matrix syntax above or the jagged matrix syntax below:

 double[][] matrix = 
 {
    new double[] { 1, 2 },
    new double[] { 3, 4 },
    new double[] { 5, 6 },
};

Special purpose matrices can also be created through specialized methods. Those include

Vector of indices

int[] idx = Matrix.Indices(0, 10);  // { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 }

Interval vector with fixed steps

double[] interval = Matrix.Interval(from: -2, to: 4); // { -2, -1, 0, 1, 2, 3, 4 };

Important matrices

double[,] I = Matrix.Identity(3);     // creates a 3x3 identity matrix
double[,] magic = Matrix.Magic(5);    // creates a magic square matrix of size 5

double[] ones = Matrix.Vector(5, 1.0);   // generates { 1, 1, 1, 1, 1 }
double[,] diagonal = Matrix.Diagonal(v); // matrix with v on its diagonal

Matrix operations

Albeit being simple double[] matrices, the framework leverages .NET extension methods to support all basic matrix operations.

double[] vector = { 0, 2, 4 };
double[] a = vector.ElementwiseMultiply(2); // vector .* 2, generates { 0,  4,  8 }
double[] b = vector.ElementwiseDivide(2);   // vector ./ 2, generates { 0,  1,  2 }
double[] c = vector.ElementwisePower(2);    // vector .^ 2, generates { 0,  4, 16 }

You can also perform operations between vectors and matrices:

  // Declare two vectors
  double[] u = { 1, 6, 3 };
  double[] v = { 9, 4, 2 };

  // Products between vectors
  double     inner     = u.InnerProduct(v);      // 39.0
  double[,]  outer     = u.OuterProduct(v);      // 
  double[]   kronecker = u.KroneckerProduct(v);  // { 9, 4, 2, 54, 24, 12, 27, 12, 6 }
  double[][] cartesian = u.CartesianProduct(v);  // all possible pair-wise combinations

  // Addition
  double[] addv = u.Add(v); // { 10, 10, 5 }
  double[] add5 = u.Add(5); // {  6, 11, 8 }

  // Elementwise operations
  double[] abs = u.Abs();   // { 1, 6, 3 }
  double[] log = u.Log();   // { 0, 1.79, 1.09 }

  // Apply *any* function to all elements in a vector
  double[] cos = u.Apply(Math.Cos); // { 0.54, 0.96, -0.989 }
  u.ApplyInPlace(Math.Cos); // can also do optionally in-place

Operations involving a matrix

  // Declare a matrix
  double[,] M = 
  {
     { 0, 5, 2 },
     { 2, 1, 5 }
  };

  // Extract a subvector from v:
  double[] vcut = v.Submatrix(0, 1); // { 9, 4 }

  // Some operations between vectors and matrices
  double[] Mv = m.Multiply(v);    //  { 24, 32 }
  double[] vM = vcut.Multiply(m); // { 8, 49, 38 }

  // Some operations between matrices
  double[,] Md = m.MultiplyByDiagonal(v);   //   { { 0, 20, 4 }, { 18, 4, 10 } }
  double[,] MMt = m.MultiplyByTranspose(m); //   { { 29, 15 },   { 15, 30 }    }

Please note this is by no means an extensive list; please take a look on all members available on this class or (preferably) use IntelliSense to navigate through all possible options when trying to perform an operation.

Decompositions

One common task in matrix manipulation is to decompose a matrix into various forms. Some examples of the decompositions supported by the framework are listed below. Those decompositions can be used to solve linear systems, compute matrix inverses and pseudo-inverses and extract other useful information about data.

Decompositions	Multidimensional	Jagged
Cholesky	(double)(float)(decimal)	(double)(float)(decimal)
Eigenvalue (EVD)	(double)(float)
Generalized Eigenvalue [1]	(double)
Nonnegative Factorization	(double)
LU	(double)(float)(decimal)	(double)(float)(decimal)
QR	(double)(float)(decimal)
Singular value (SVD)	(double)(float)

Loading and parsing from text

The framework can load any MATLAB-compatible .mat file using the MatReader class. It can also parse matrices written in MATLAB, Octave, Mathematica and C#/NET formats directly from text. Examples are shown below.

// From free-text format
double[,] a = Matrix.Parse(@"1 2
                             3 4");

// From MATLAB/Octave format
matrix[,] b = Matrix.Parse("[1 2; 3 4]", OctaveMatrixFormatProvider.InvariantCulture);

// From C# multi-dimensional array format
string str = @"double[,] matrix = 
               {
                  { 1, 2 },
                  { 3, 4 },
                  { 5, 6 },
               }";

double[,] c = Matrix.Parse(str, CSharpMatrixFormatProvider.InvariantCulture);

Optimization

Another very common task in mathematics, machine learning and statistics is the need to optimize (i.e. find the maximum or minimum) of a function. The Accord.Math.Optimization namespace contains classes for constrained and unconstrained optimization, such as for example (Conjugate Gradient) (CG), Broyden–Fletcher–Goldfarb–Shanno (BFGS) and the Bounded Broyden–Fletcher–Goldfarb–Shanno (Bounded-BFGS) for unconstrained optimization; the Goldfarb-Idnani solver for Quadratic Programming (QP) problems, the Augmented Lagrangian method for nonlinear optimization and gradient-free methods such as Cobyla. It also provides methods for univariate search such as Brent's method for finding function roots, maximum or minimum values of an univariate function.

To find the root of an univariate function, you can use

Func<double, double> function = x => x * x * x + 2 * x * x - 10 * x;

// And now we can create the search algorithm:
BrentSearch search = new BrentSearch(function, -4, 3);

// Finally, we can query the information we need
double max = search.Maximize();  // occurs at -2.61
double min = search.Minimize();  // occurs at  1.27
double root = search.FindRoot(); // occurs at  0.50

For multivariate functions, you use

// Suppose we would like to find the minimum of the function
// 
//   f(x,y)  =  -exp{-(x-1)²} - exp{-(y-2)²/2}
// 

// First we need write down the function either as a named
// method, an anonymous method or as a lambda function:

Func<double[], double> f = (x) =>
    -Math.Exp(-Math.Pow(x[0] - 1, 2)) - Math.Exp(-0.5 * Math.Pow(x[1] - 2, 2));

// Now, we need to write its gradient, which is just the
// vector of first partial derivatives del_f / del_x, as:
// 
//   g(x,y)  =  { del f / del x, del f / del y }
// 

Func<double[], double[]> g = (x) => new double[] 
{
    // df/dx = {-2 e^(-    (x-1)^2) (x-1)}
    2 * Math.Exp(-Math.Pow(x[0] - 1, 2)) * (x[0] - 1),

    // df/dy = {-  e^(-1/2 (y-2)^2) (y-2)}
    Math.Exp(-0.5 * Math.Pow(x[1] - 2, 2)) * (x[1] - 2)
};

// Finally, we can create the L-BFGS solver, passing the functions as arguments
var lbfgs = new BroydenFletcherGoldfarbShanno(numberOfVariables: 2, function: f, gradient: g);

// And then minimize the function:
bool success = lbfgs.Minimize();
double minValue = lbfgs.Value;
double[] solution = lbfgs.Solution;

If you have a quadratic programming problem with constraints, you can use

// Solve the following optimization problem:
// 
//  max f(x) = -2x² + xy - y² + 5y
// 
//  s.t.   x - y  ==   5  (x minus y should be equal to 5)
//             x  >=  10  (x should be greater than or equal to 10)
// 
// 

// Create our objective function using a text string
var f = new QuadraticObjectiveFunction("-2x² + xy - y² + 5y");

// Now, create the constraints
List<LinearConstraint> constraints = new List<LinearConstraint>();
constraints.Add(new LinearConstraint(f, "x - y ==  5"));
constraints.Add(new LinearConstraint(f, "    x >= 10"));

// Now we create the quadratic programming solver for 2 variables, using the constraints.
GoldfarbIdnani solver = new GoldfarbIdnani(f, constraints);

// And attempt to solve it.
double maxValue = solver.Maximize();

If you have a general optimization problem with constraints, then you can use

// Suppose we would like to minimize the following function:
// 
//    f(x,y) = min 100(y-x²)²+(1-x)²
// 
// Subject to the constraints
// 
//    x >= 0  (x must be positive)
//    y >= 0  (y must be positive)
// 

// In this example we will be using some symbolic processing. 
// The following variables could be initialized to any value.

double x = 0, y = 0;


// First, we create our objective function
var f = new NonlinearObjectiveFunction(

    // This is the objective function:  f(x,y) = min 100(y-x²)²+(1-x)²
    function: () => 100 * Math.Pow(y - x * x, 2) + Math.Pow(1 - x, 2),

    // The gradient vector:
    gradient: () => new[] 
    {
        2 * (200 * Math.Pow(x, 3) - 200 * x * y + x - 1), // df/dx = 2(200x³-200xy+x-1)
        200 * (y - x*x)                                   // df/dy = 200(y-x²)
    }

);


// Now we can start stating the constraints
var constraints = new List<NonlinearConstraint>();

// Add the non-negativity constraint for x
constraints.Add(new NonlinearConstraint(f,

    // 1st constraint: x should be greater than or equal to 0
    function: () => x, shouldBe: ConstraintType.GreaterThanOrEqualTo, value: 0,

    gradient: () => new[] { 1.0, 0.0 }
));

// Add the non-negativity constraint for y
constraints.Add(new NonlinearConstraint(f,

    // 2nd constraint: y should be greater than or equal to 0
    function: () => y, shouldBe: ConstraintType.GreaterThanOrEqualTo, value: 0,

    gradient: () => new[] { 0.0, 1.0 }
));


// Finally, we create the non-linear programming solver
var solver = new AugmentedLagrangianSolver(2, constraints);

// And attempt to solve the problem
double minValue = solver.Minimize(f);

And finally, if you have a gradient that is difficult to compute or hard to derive analytically, you can use a gradient-free optimization algorithm with support for optional constraints as

// We would like to find the minimum of min f(x) = 10 * (x+1)^2 + y^2
Func<double[], double> function = x => 10 * Math.Pow(x[0] + 1, 2) + Math.Pow(x[1], 2);
   
// Create a cobyla method for 2 variables
Cobyla cobyla = new Cobyla(2, function);

bool success = cobyla.Minimize();

double minimum = minimum = cobyla.Value; // Minimum should be 0.
double[] solution = cobyla.Solution;     // Vector should be (-1, 0)

or in the case you have a constrained problem, you can use

// We will optimize the 2-variable function f(x, y) = -x -y
var f = new NonlinearObjectiveFunction(2, x => -x[0] - x[1]);

// Under the following constraints
var constraints = new[]
{
   new NonlinearConstraint(2, x =>             x[1] - x[0] * x[0] >= 0),
   new NonlinearConstraint(2, x =>  1 - x[0] * x[0] - x[1] * x[1] >= 0),
};

// Create a Cobyla algorithm for the problem
var cobyla = new Cobyla(function, constraints);

// Optimize it
bool success = cobyla.Minimize();
double minimum = cobyla.Value;        // Minimum should be -2 * sqrt(0.5)
double[] solution = cobyla.Solution;  // Vector should be [sqrt(0.5), sqrt(0.5)]

Derivatives and Integrals

Special functions

Complete sets of special functions, such as factorials, log-factorials, combinatorics and more specialized functions, such as the [http://accord.googlecode.com/svn/docs/html/AllMembers_T_Accord_Math_Bessel.htm Bessel], [http://accord.googlecode.com/svn/docs/html/AllMembers_T_Accord_Math_Gamma.htm Gamma], [http://accord.googlecode.com/svn/docs/html/AllMembers_T_Accord_Math_Beta.htm Beta] and [http://accord.googlecode.com/svn/docs/html/AllMembers_T_Accord_Math_Normal.htm Normal (Gaussian)] functions.