F# INumerics Interface, and Matrix class

The blog post is inspired by a question on StackOverflow. It asks whether the Matrix class in F# PowerPack supports a user defined type.

First the answer is: YES.

But supporting that is non-trivial. I found very few material online on this subject, so I decided to to write a blog post on it. Since matrix is the very core concept in data mining and machine learning, our digging into the F# Matrix and relevant interfaces would also be valuable.

In this post, I will lead to read a part of the F# PowerPack source code, which is in source\fsppack\FSharp.PowerPack\math folder.

Although the latest version is 1.9.9.9, I would suggest to use the source code of 1.9.7.8, which contains a lapack subfolder in the math folder. What’s this? The researchers in MSR did a great wrapper to enable F# to call common linear algebra libraries, e.g. the free lapack and the commercial IMK (Intel Math Kernel Library). However, this part of code is removed from the 1.9.9.9 release. These code will be merged into Math.Net project, which is still under development. This is just a kind notice, we don’t need to know lapack wrapper today.

Let’s get back to that question. In F#, we can use matrices similar to Matlab:

let B = matrix [ [ 1.0; 7.0 ];
               [ 1.0; 3.0 ] ]

where matrix is a actually a function, its declaration/signature is

val matrix : seq<#seq<float>> –> matrix

and the implementation is:

let matrix ll = Matrix.ofSeq ll

As a list implements the IENumerable interface, so we can use a list to initialize a matrix, so is array.

Let see the type of B:

val B : matrix = matrix [[1.0; 7.0]
                       [1.0; 3.0]]

So we can see that matrix is also a type name, it is ok to have the same name for a type and a function. In the source file Matrix.fs, we can find the definition of type matrix:

type matrix = Matrix<float>

It is defined as Matrix<float>, where Matrix<_> is called a meta type or polymorphic type. Meta types are used to create other type by putting a element type in between <>. You can think it as an equivalent to C++ templates, or Java/.Net Generics. (There are differences, however that would be a long story.)

Also notice that Matrix in Matrix<_> is a meta type, while Matrix in Matrix.ofSeq is a module name. So they also can have the same name as they are different things.

One thing we can come to mind is that, we can define our integer matrix type as

type intmatrix = Matrix<int>

and we also need a constructor for it

let intmatrix ll = Matrix.Generic.ofSeq  ll

Matrix.Generic module contains functions for generic type matrices, while the functions in Matrix module are only for float matrices.

We can happily do some operations on integer matrices:

let C = intmatrix [ [1;2]; [3;4] ];
let D = intmatrix [ [1;1]; [1;1] ];
let E = Matrix.Generic.inplaceAdd C D
C+D
C * 10

So far so good. So we can put float, int, etc into Matrix<_>, can we put:

type Pixel = 
  | P of int*int*int

into it? Kind of yes:

let C = pmatrix [ [P(2,2,2);P(2,2,2)]; [P(2,2,2);P(2,2,2)] ];

But when adding two such matrices, of course error occurs!

We must lack something, e.g. at least we need to define the behavior of the add operator for Pixel type.

The idea is that every type that could put into Matrix<_> is associated with an instance of INumeric interface. Notice that I use the word ‘associate’. It means that for each number type T, we associate an object with it, which tells the matrix class how to add, subtract, multiply, etc. on T. E.g. we have Int32Numerics defined in INumeric.fs, which is associated with int type:

let Int32Numerics =
{ new IIntegral<int32> with
     member __.Zero = 0
     member __.One = 1
     member __.Add(a,b) = a + b
     member __.Subtract(a,b) = a - b
     member __.Multiply(a,b) = a * b
     member __.Equals(a,b) = (a = b)
     member __.Compare(a,b) = compare a b
     member __.Negate(a) = - a
     member __.Abs(a) = a
     member __.ToBigInt(a) = new BigInteger(a)
     member __.OfBigInt(a) = int32 a
     member __.Sign(a) = Math.Sign(a)
     member __.Modulus(a,b) = a % b
     member __.Divide(a,b) = a / b
     member __.DivRem(a,b) = (a / b, a % b)
     member __.ToString((x:int32),fmt,fmtprovider) =
            x.ToString(fmt,fmtprovider)
     member __.Parse(s,numstyle,fmtprovider) =
            System.Int32.Parse(s,numstyle,fmtprovider)
  interface INormFloat<int32> with
     member __.Norm(x) = float (abs x)
}

where IIntegral<> is an interface inherited from INumeric<>. 

Let’s move on to user defined types. We know there is a Big Rational type defined in F# PowerPack. We can easily use it as:

let rmatrix ll = Matrix.Generic.ofSeq  ll
let C = rmatrix [ [1N;2N]; [3N;4N] ];
let D = rmatrix [ [1N;1N]; [1N;1N] ];
C+D

The rational number is defined in q.fs. However, the code to enable it to be used in Matrix<_> are in other places. First, In INumerics.fs, we find:

let BigNumNumerics =
{ new IFractional<bignum> with
     member __.Zero = BigNum.Zero
     member __.One = BigNum.One
     member __.Add(a,b)      = a + b
     member __.Subtract(a,b) = a - b
     member __.Multiply(a,b) = a * b
     member __.Equals(a,b) = (a = b)

How will Matrix know a number type is associated with its numeric? The code is in associations.fs:

let ht =
  let ht = new System.Collections.Generic.Dictionary<Type,obj>()
  let optab =
      [ typeof<float>,   (Some(FloatNumerics    :> INumeric<float>) :> obj);
        typeof<int32>,   (Some(Int32Numerics    :> INumeric<int32>) :> obj);
        typeof<int64>,   (Some(Int64Numerics    :> INumeric<int64>) :> obj);
        typeof<BigInteger>,  (Some(BigIntNumerics   :> INumeric<BigInteger>) :> obj);
        typeof<float32>, (Some(Float32Numerics  :> INumeric<float32>) :> obj);
        typeof<Microsoft.FSharp.Math.Complex>, (Some(ComplexNumerics :> INumeric<Microsoft.FSharp.Math.Complex>) :> obj);
        typeof<bignum>,  (Some(BigNumNumerics   :> INumeric<bignum>) :> obj); ]
   
  List.iter (fun (ty,ops) -> ht.Add(ty,ops)) optab;
  ht

The types above are automatically supported. If we define a new type, we need to add the association of the type and its numeric into the variable ht.

By now, we know all the secrets. We can code our Pixel type now:

type Pixel =

    | P of int*int*int
  static member makePixel (a,b,c) = P (a,b,c)

  override p.ToString() =
      let (P(x,y,z)) = p
      "("+x.ToString()+", "+y.ToString()+", "+z.ToString()+")"
    
  static member (+) (P(a,b,c), P(x,y,z)) = P (a+x, b+y, c+z)

  static member (-) (P(a,b,c), P(x,y,z)) = P (a-x, b-y, c-z)
  static member (*) (P(a,b,c), P(x,y,z)) = P (0, 0, 0)
  static member (/) (P(a,b,c), P(x,y,z)) = P (0, 0, 0)
  static member (~+) (P(a,b,c), P(x,y,z)) = P (0, 0, 0)
  static member (~-) (P(a,b,c), P(x,y,z)) = P (0, 0, 0)

Its numerics:

let PixelNumerics =
  { new INumeric<Pixel> with
      member op.Zero = P(0,0,0)
      member op.One = P(1,1,1)
      member op.Add(a,b) = a + b
      member op.Subtract(a,b) = a - b
      member op.Multiply(a,b) = a * b
      member ops.Negate(a) = a
      member ops.Abs(a) = a
      member ops.Equals(a, b) = false
      member ops.Compare(a, b) = 0
      member ops.Sign(a) = failwith "not defined."
      member ops.ToString(x,fmt,fmtprovider) =  failwith "not implemented"
      member ops.Parse(s,numstyle,fmtprovider) = failwith "not implemented"
  }

add the association:

GlobalAssociations.RegisterNumericAssociation(PixelNumerics)