NLinear: Generic linear algebra library for .NET

Publié dans: 

NLinear is a generic linear algebra library for .NET. Probably It’s not the most complete nor the most efficient open source .NET linear algebra toolkit. But the library has a distinguishing feature: genericity. In fact NLinear isn’t bound to a specific scalar type (for example matrix of float or of double) instead it can operate on any type satisfying some requirements. 

Back to basics: Vector spaces

Linear algebra deals with vector spaces.A familiar example of a vector space is the n-dimensional euclidean space. A vector space is stable under vectors addition and scalar multiplication. Addition and multiplication must satisfy certain conditions.

Vector spaces are defined over a scalar field. A field is a particular type of rings, an algebraic structure defined as :

A ring is a set A with two laws of composition called respectively addition and multiplication, satisfying the following axioms:

  1. Under addition A is a commutative group.
  2. Multiplication is associative and possess an identity element.
  3. Multiplication is distributive with respect to addition.

A ring K is called a field if it does not consist only of 0 and every non-zero element of K is invertible.

The objective is to create a library that works with arbitrary scalar fields. In the next section we try to define the scalar field in C#.

Modeling a field

The syntax of the library is based on C# generics. Matrices and Vectors have a type parameter. To create a complex matrix it suffice to write:

Matrix<Complex> m = new Matrix<Complex>();

For conventional real matrix we write:

Matrix<double> m = new Matrix<double>();

The first solution I come up with was to use constraints on type parameter, to ensure that the type implements an interface representing a Field:

class Matrix<T> where T : IField
{
   (...)
}

    interface IField
    {

        // The first law of composition (addition)
        IField Plus(IField x);

       
        // The second law of composition (multiplication)
        IField Times(IField x);

       
        // Additive identity (0 + x = x + 0 = x)
        IField AdditiveIdentity();

       
        // Multiplicative identity (x * 1 = 1 * x = x)
        IField MultiplicativeIdentity();

       (...)

    }

This implementation have two fundamental shortcomings.To begin with, we can’t define operators in C# interfaces. Thus if a, b are instance of a type implementing IField, we can’t write a + b. Instead we should write a.Plus(b). And the other big limitation is that the implementation doesn’t work for basic numeric types such as int, double, float… because they doesn’t implement IField. Another ways must be explored.

There is some fundamental differences between C++ templates and C# generics, Probably It would be futile to try to implement a solution inspired by a generic C++ implementation. Instead I looked for ways to circumvent C# generics limitations. I found a particular usage of C# code expression tree to implement generic operators discovered by Marc Gravell[1].

Expression tree is a feature introduced in C# 3.0 mostly used by LINQ and the DLR. It enables the representation of code in a tree structure. One could build, modify, compile and run the code represented as a tree in run-time.

Given the expression (1+1) it could be represented as an Add expression:

Expression sumExpr = Expression.Add(
      Expression.Constant(1),
      Expression.Constant(1)
);

The above expression could be compiled to a lambda:

var compiledExper = Expression.Lambda>(sumExpr).Compile();

If we execute the compiled expression we obtain the result 2.

Console.WriteLine(compiledExper());

We now have all the ingredients to define the class Numeric<T>. This class represent a type with a law of composition (+).

    class Numeric<T>
    {
        T data;

        static Func<T, T, T> add = null;    

        public Numeric(T data)
        {

            this.data = data;

            if (add == null)
            {
                ParameterExpression argX = Expression.Parameter(typeof(T), "x");
                ParameterExpression argY = Expression.Parameter(typeof(T), "y");
              
                BinaryExpression body = Expression.Add(argX, argY);

                add = Expression.Lambda<Func<T, T, T>>(body, argX, argY).Compile();
            }
        }

       public static Numeric<T> operator+ (Numeric<T> x, Numeric<T> y)
        {
            T result = add(x.GetValue(), y.GetValue());
            return new Numeric<T>(result);
        }

        public T GetValue()
        {
            return data;
        }

        public override string ToString()
        {
            return string.Format("{0}", data);
        }
    }

To perform addition between two integers we can use :

   Numeric<int> a = new Numeric<int>(1);
            Numeric<int> b = new Numeric<int>(1);

            Console.WriteLine(a + b);

Suppose that we pass a type that doesn’t have addition as a law of composition:

  Numeric<DateTime> a = new Numeric<DateTime>(new DateTime(0));
            Numeric<DateTime> b = new Numeric<DateTime>(new DateTime(8));

What happens when the calculation a + b is invoked ?

Console.WriteLine(a + b);

The above code compiles ! but when executing it an exception of the type InvalidOperationException is thrown. The content of this exception: “The binary operator Add is not defined for the type DateTime”.

Numeric<T> could be completed to include operators *, / et - and so on. A very handy tool is the Operator class defined in the library MiscUtil[2] that provides an implementation of all the operators.

Let’s define a simplified version of a 2d vector. Instead of working directly on a numeric type T, NLinear uses Numeric<T> to access to it’s operators. We also define an implicit cast from T to Numeric<T> thus we could use the unmodified algorithm operating on type T.

 /// <summary>
    /// A 2-dimensional vector/point
    /// </summary>
    /// <typeparam name="T">A numeric type (float, int, ...)</typeparam>
    public struct Vec2<T> : IEquatable<Vec2<T>>
        where T : IEquatable<T>
    {
        /// <summary>
        ///  The vector components.
        /// </summary>
        Numeric<T> x, y;
       
        /// <summary>
        /// Constructor (x y) = (a b)
        /// </summary>
        /// <param name="a"></param>
        public Vec2(T a, T b)
        {
            x = a;
            y = b;
        }

        /// <summary>
        /// The dot product.
        /// </summary>
        /// <param name="v"></param>
        /// <returns></returns>
        public T Dot(Vec2<T> v)
        {
            return x * v.x + y * v.y;
        }

        /// <summary>
        /// Cross-product
        /// </summary>
        /// <param name="v"></param>
        /// <returns></returns>
        public T Cross(Vec2<T> v)
        {
            return x * v.y - y * v.x;
        }
    }

Using the above defined class is straightforward :

         //Declare two unit vectors e1, e2
           Vec2<BigInteger> e1 = new Vec2<BigInteger>(1, 0);
           Vec2<BigInteger> e2 = new Vec2<BigInteger>(0, 1);

           Vec2<BigInteger> e3 = e1.Cross(e2);

Using this approach for all structures having a type parameter leads to generic implementation. Code expression is compiled only at once thus there is virtually no major performance penalties. But the implementation has also some limitations.

We are constantly looking for ways to enhance NLinear, you can contribute code or ideas on http://nlinear.codeplex.com/ and IP-TECH Labs, or just by typing your comment below.


[1]http://www.yoda.arachsys.com/csharp/genericoperators.html

[2]http://www.yoda.arachsys.com/csharp/miscutil/