Mathematics
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.
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
int[] idx = Matrix.Indices(0, 10); // { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 }double[] interval = Matrix.Interval(from: -2, to: 4); // { -2, -1, 0, 1, 2, 3, 4 };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// 1.1 Using standard .NET declaration double[,] A = { {1, 2, 3}, {6, 2, 0}, {0, 0, 1} };
double[,] B = { {2, 0, 0}, {0, 2, 0}, {0, 0, 2} };
{ // 1.2 Using Accord extension methods double[,] Bi = Matrix.Identity(3).Multiply(2); double[,] Bj = Matrix.Diagonal(3, 2.0); // both are equal to B
// 1.2 Using Accord extension methods with implicit typing
var I = Matrix.Identity(3);
}
Albeit being simple double[] matrices, the framework leverages .NET extension methods to support all basic matrix operations.
double[] v = { 4, 5, 6 };
double[] a = v.ElementwiseMultiply(2); // v .* 2: { 8, 10, 12 }
double[] b = v.ElementwiseDivide(2); // v ./ 2: { 2, 2.5, 3 }
double[] c = v.ElementwisePower(2); // v .^ 2: { 16, 25, 36 }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-placeOperations 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.
// 3.1 Calculating the determinant
double det = A.Determinant();
// 3.2 Calculating the trace
double tr = A.Trace();
// 3.3 Computing the sum vector
{
double[] sumVector = A.Sum();
// 3.3.1 Computing the total sum of elements
double sum = sumVector.Sum();
// 3.3.2 Computing the sum along the rows
sumVector = A.Sum(0); // Equivalent to Octave's sum(A, 1)
// 3.3.2 Computing the sum along the columns
sumVector = A.Sum(1); // Equivalent to Octave's sum(A, 2)
}{
// 2.1 Addition
var C = A.Add(B);
// 2.2 Subtraction
var D = A.Subtract(B);
// 2.3 Multiplication
{
// 2.3.1 By a scalar
var halfM = A.Multiply(0.5);
// 2.3.2 By a vector
double[] m = A.Multiply(new double[] { 1, 2, 3 });
// 2.3.3 By a matrix
var M = A.Multiply(B);
// 2.4 Transposing
var At = A.Transpose();
}
}
// 2.5 Elementwise operations
// 2.5.1 Elementwise multiplication
A.ElementwiseMultiply(B); // A.*B
// 2.5.1 Elementwise division
A.ElementwiseDivide(B); // A./B// 4.1 Computing the inverse
var invA = A.Inverse();
// 4.2 Computing the pseudo-inverse
var pinvA = A.PseudoInverse();
// 4.3 Solving a linear system (Ax = B)
var x = A.Solve(B);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) |
Below are some examples on how to use the decompositions.
// Singular value decomposition
{
SingularValueDecomposition svd = new SingularValueDecomposition(A);
var U = svd.LeftSingularVectors;
var S = svd.Diagonal;
var V = svd.RightSingularVectors;
}
// or (please see documentation for details)
{
SingularValueDecomposition svd = new SingularValueDecomposition(A.Transpose());
var U = svd.RightSingularVectors;
var S = svd.Diagonal;
var V = svd.LeftSingularVectors;
}
// Eigenvalue decomposition
{
EigenvalueDecomposition eig = new EigenvalueDecomposition(A);
var V = eig.Eigenvectors;
var D = eig.DiagonalMatrix;
}
// QR decomposition
{
QrDecomposition qr = new QrDecomposition(A);
var Q = qr.OrthogonalFactor;
var R = qr.UpperTriangularFactor;
}
// Cholesky decomposition
{
CholeskyDecomposition chol = new CholeskyDecomposition(A);
var R = chol.LeftTriangularFactor;
}
// LU decomposition
{
LuDecomposition lu = new LuDecomposition(A);
var L = lu.LowerTriangularFactor;
var U = lu.UpperTriangularFactor;
}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);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.50For 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)]You can compute your own derivatives / gradient vectors using numerical methods, such as Newton's Finite Differences method.
// Create a simple function with two parameters: f(x,y) = x² + y
Func<double[], double> function = x => Math.Pow(x[0], 2) + x[1];
// Create a new finite differences calculator
var calculator = new FiniteDifferences(2, function);
// The gradient function should be g(x,y) = <2x, 1>
Func<double[], double[]> gradient = calculator.Compute;
// Evaluate the gradient function at the point (2, -1)
double[] result = gradient(2, -1); // answer is (4, 1)For integration methods, there are several approaches you can use. You can use the simplistic but efficient Trapezoid approximation:
or more elaborated methods such as InfiniteAdaptiveGaussKronrod or even Monte-Carlo approximation methods
Check the other methods available in the Accord.Math.Differentiation and Accord.Math.Integration namespaces.
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.
#region 1. Declaring matrices
// 1.1 Using standard .NET declaration
double[,] A =
{
{1, 2, 3},
{6, 2, 0},
{0, 0, 1}
};
double[,] B =
{
{2, 0, 0},
{0, 2, 0},
{0, 0, 2}
};
{
// 1.2 Using Accord extension methods
double[,] Bi = Matrix.Identity(3).Multiply(2);
double[,] Bj = Matrix.Diagonal(3, 2.0); // both are equal to B
// 1.2 Using Accord extension methods with implicit typing
var I = Matrix.Identity(3);
}
#endregion
#region 2. Matrix Operations
{
// 2.1 Addition
var C = A.Add(B);
// 2.2 Subtraction
var D = A.Subtract(B);
// 2.3 Multiplication
{
// 2.3.1 By a scalar
var halfM = A.Multiply(0.5);
// 2.3.2 By a vector
double[] m = A.Multiply(new double[] { 1, 2, 3 });
// 2.3.3 By a matrix
var M = A.Multiply(B);
// 2.4 Transposing
var At = A.Transpose();
}
}
// 2.5 Elementwise operations
// 2.5.1 Elementwise multiplication
A.ElementwiseMultiply(B); // A.*B
// 2.5.1 Elementwise division
A.ElementwiseDivide(B); // A./B
#endregion
#region 3. Matrix characteristics
{
// 3.1 Calculating the determinant
double det = A.Determinant();
// 3.2 Calculating the trace
double tr = A.Trace();
// 3.3 Computing the sum vector
{
double[] sumVector = A.Sum();
// 3.3.1 Computing the total sum of elements
double sum = sumVector.Sum();
// 3.3.2 Computing the sum along the rows
sumVector = A.Sum(0); // Equivalent to Octave's sum(A, 1)
// 3.3.2 Computing the sum along the columns
sumVector = A.Sum(1); // Equivalent to Octave's sum(A, 2)
}
}
#endregion
#region 4. Linear Algebra
{
// 4.1 Computing the inverse
var invA = A.Inverse();
// 4.2 Computing the pseudo-inverse
var pinvA = A.PseudoInverse();
// 4.3 Solving a linear system (Ax = B)
var x = A.Solve(B);
}
#endregion
#region 5. Special operators
{
// 5.1 Finding the indices of elements
double[] v = { 5, 2, 2, 7, 1, 0 };
int[] idx = v.Find(e => e > 2); // finding the index of every element in v higher than 2.
// 5.2 Selecting elements by index
double[] u = v.Submatrix(idx); // u is { 5, 7 }
// 5.3 Converting between different matrix representations
double[][] jaggedA = A.ToArray(); // from multidimensional to jagged array
// 5.4 Extracting a column or row from the matrix
double[] a = A.GetColumn(0); // retrieves the first column
double[] b = B.GetRow(1); // retrieves the second row
// 5.5 Taking the absolute of a matrix
var absA = A.Abs();
// 5.6 Applying some function to every element
var newv = v.Apply(e => e + 1);
}
#endregion
#region 7. Vector operations
{
double[] u = { 1, 2, 3 };
double[] v = { 4, 5, 6 };
var w1 = u.InnerProduct(v);
var w2 = u.OuterProduct(v);
var w3 = u.CartesianProduct(v);
double[] m = { 1, 2, 3, 4 };
double[,] M = Matrix.Reshape(m, 2, 2);
}
#endregion
#region Decompositions
{
// Singular value decomposition
{
SingularValueDecomposition svd = new SingularValueDecomposition(A);
var U = svd.LeftSingularVectors;
var S = svd.Diagonal;
var V = svd.RightSingularVectors;
}
// or (please see documentation for details)
{
SingularValueDecomposition svd = new SingularValueDecomposition(A.Transpose());
var U = svd.RightSingularVectors;
var S = svd.Diagonal;
var V = svd.LeftSingularVectors;
}
// Eigenvalue decomposition
{
EigenvalueDecomposition eig = new EigenvalueDecomposition(A);
var V = eig.Eigenvectors;
var D = eig.DiagonalMatrix;
}
// QR decomposition
{
QrDecomposition qr = new QrDecomposition(A);
var Q = qr.OrthogonalFactor;
var R = qr.UpperTriangularFactor;
}
// Cholesky decomposition
{
CholeskyDecomposition chol = new CholeskyDecomposition(A);
var R = chol.LeftTriangularFactor;
}
// LU decomposition
{
LuDecomposition lu = new LuDecomposition(A);
var L = lu.LowerTriangularFactor;
var U = lu.UpperTriangularFactor;
}
}
#endregion