F# and Data Mining

Any numerical computing environment on Java platform?

2011-10-31T19:41:00.001+08:00

Recently Tomas Petricek and I have co-written a book chapter, Numerical Computing in F#. It is a free online book chapter for his Manning book: Real World Functional Programming. This book chapter is a survey on numerical computing for .Net languages, specially for F#. I will write another article on this book chapter and present some materials that were cut down from final version on MSDN.

But today, I am looking for things on our competitor’s side -- numerical computing environments for Java platform. I am just wondering what Java guys are doing for numerical computing.

Looking at the mature numerical systems such as Matlab, Mathematica and S-plus/R, we can summarize three general elements of such a system:

1. Math library. Math library is a must for serious numerical computing. It is hard to build well-tested numerical procedures. This kind of low level stuff, e.g. linear algebra routines and FFT, is usually available as a reusable software component.

2. Interactive/Console. The Read-eval-print loop (REPL) is very important for exploratory analysis because it gives a quick feedback on some small fragments of code.

3. Plot/chart library. A plot explains what’s going on clearly and quickly!

I did some web search and found quite a few blog articles; I want to mention articles by Mikio L. Braun. He also maintains his jblas library, which is a JNI wrapper to ATLAS, the highly-optimized BLAS library.

I am aware Incanter, which is a Clojure project aiming to be R in Java platform. But I don't think rewriting basic things, e.g. PCA and linear regression, for Java Platform is a good idea since a stable implementation for these algorithms usually take years because there is just too much details in numerical algorithms.

This time, I found a new library, ScalaLab, which is a on-going project for numerical computing in Scala.

Let’s first talk the library a little bit. For matrix and linear algebra, there are also two quite stable Java libraries: colt and JAMA. The linear algebra functionality of the two libraries should be intact. But the performance is quite bad compared to native code, even several times slower to pure .Net code. [The comparison to .Net code is my personal experience. ] But one should not care too much about performance as the bottleneck of a numerical project differs from projects to projects. A working system is the most important thing, various methods could help to improve the performance. Btw, the performance difference is usually not noticeable when we work on small or sample datasets. I use R quite often recently, the matrix library associated with its Windows version is from Netlib, a normal native performance one. I didn’t have any performance complain over months.

For REPL, basically it is a language issue. After reading through online, I think I will like something built on Scala, such as the ScalaLab environment. Among the mature/quite-mature JVM languages, Scala is most similar to F#. The advantage of Scala to F# is that Scala is a component language, which means Scala has better OO features for library designers and software architects. However, on the small/low level of programming (e.g. a small function), F# feels much more functional and pure than Scala. This, I think, is partially because F# “IS” an ML, while Scala borrows some syntax from Java and C#. The other reason is that OO is at the heart of Scala – everything is an object. Scala’s OO model is greatly influenced by Smalltalk’s. I haven’t learned any Smalltalk. But after some learning in Scala, I can appreciate why Alan Kay, the inventor of Smalltalk, says, “Actually I made up the term "object-oriented", and I can tell you I did not have C++ in mind. “. So even functional programming in Scala is OO-emulated, while F# has a functional base first and then the implementation using .Net is transparent to its syntax. Ok… the above comparison is just off the topic. Anyway, Scala can serve as a good static script language with a REPL.

For the plotting library, JFreeChart seems to be the standard. It is also the plotting library used in Incanter [For this part, I think Incanter really does a great job!]. I am not sure whether there are some advanced or commercial plotting libraries for Java. But I know that .Net has tons of that! Anyway, JFreeChart is ok for daily use.

The last thing, which I haven’t listed above, is parallel computing! Java platform, or more specific, the Java Virtual Machine, proves to be one of the most stable multithreaded systems ever built. For example, Hadoop is written in Java. Although people in Google said that Hadoop is way slow compared to their C++ MapReduce. A lot of big companies are using it. Two months ago, I read the following nice paper by Martin Odersky and Kunle Olukotun:

Language Virtualization for Heterogeneous Parallel Computing

It is just published and already has 28 citations on Google Scholar. The idea is to build a DSL, which serves a middle ware, to write high level programs, e.g. equations involving some matrix operations. Underlying the DSL, there are complicated optimizations going on to transform it and make it parallel using heterogeneous technologies, e.g. GPU and multi-core. So the programmer does not need to know anything about parallel computing and his program is executed in parallel utilizing the multi-cores! The DSL is an internal one built in Scala, another example showing the great expressiveness of Scala. Since such idea is not new, the novel part of the paper is to demonstrate how Scala perfectly fits the requirements of such a paradigm.

I don’t quite have a conclusion for this post. Just some random thoughts on numerical computing in Java and an advertisement for my articles on MSDN

A functional red black tree with dynamic order statistics in F#

2011-09-14T23:04:00.001+08:00

Yin Zhu

14 September 2011

(The code for this article is at codeplex.)

Binary search trees play an extremely important role in computer science. They were hot research topics back to 1970s and early 1980s. While the plain binary search tree is of little usage, balanced search trees are now being used everywhere. For example, most C++ STL implementations use Red Black Tree for its set and map containers. The F# Core library uses a functional version of AVL tree for Set and Map, which are both persistent/functional containers.

Special binary search trees also have other functionalities other than implementing standard set or map containers. For example, splay tree has an efficient split operation – splitting a dynamic set into two in O(log n) time: one is less than the given key; the other is greater than or equal to the given key. This operation could be used to optimize network flow algorithms. Details could be found in Tarjan’s classical (also tiny) data structures book published in 1983.

In this post, I’d like to implement the classical Red Black Tree in F#. I searched on the internet and found that there were some Red Black Tree code snippets in F#, but none of them are complete. Some only implemented the insert operation, which is relatively easy. Some are incorrect.

My main purpose for implementing the Red Black Tree in F# is to test whether I am still good at implementing complex data structures (I implemented quite a few classical data structures and algorithms during my early undergraduates for ACM/ICPC competitions). But there are new things: 1) the implementation language is F# rather than C++, 2) the implementation should be functional rather than imperative, and 3) the interface design should be general, and allows code reuse.

One function I always miss from standard set containers is dynamic order statistics – finding the ith smallest element in a dynamic set. The query time should be in O(log n). To accomplish this, we need add a size field to the internal tree node. Adding extra fields into a data structure is called augmenting data structures (Chapter 14, CLRS). Augmenting data structures is useful when a function is not built in a classical data structure but can be added by genuinely modifying the classical data structure. It is a very important skill when solving algorithmic competition problems.

I divided this article into several sections. I will first introduce the basics of red black trees, and then we move the insertion part, and order statistics, and then the hardest part – deletion. I have a section describing the testing procedures that ensure the correctness of the implementation. At last, I will leave the specifics of Red Black Trees, and cover the tree traversal in general and show some more functional stuff: continuation and a continuation Monad.

Red Black Tree Basics

A Red Black Tree is a binary search tree with the following extra properties:

1. The root is black and the leaf nodes are black (note: leaf nodes do not contain keys, only internal nodes contain keys).

2. A red node’s two children must be black.

3. For any subtree rooted at x, all the simple paths down to the leaf nodes contain the same number of black nodes.

With the above properties, the Red Black tree is able to perform queries such as search, insert and delete all in O(log n) time.

The F# definition of the tree is given below:

type Color =
   Red | Black | DoubleBlack
  
type RBTree<'a when 'a:comparison> =
   | Node of Color * 'a * int * RBTree<'a> * RBTree<'a>
   | Leaf of Color

Notice that a node could be Red or Black. For black, there are two cases Black and DoubleBlack. The DoubleBlack color is used only when we are dealing with deletion, and in a valid Red Black tree, there would be of no DoubleBlack node. I will explain DoubleBlack in the delete section.

A Red Black Tree has two kinds of nodes: internal nodes storing keys and dummy leaf nodes. The internal node has a color field, the key filed, the size of the subtree and left and right children. The leaf node has two possible colors: Black or DoubleBlack.

Let’s write two simple recursive functions working on the tree structure to get us warmed up:

let rec contains x = function
   | Leaf _ -> false
   | Node (_, y, _, a, b) ->
       if x = y then true
       elif x < y then contains x a
       else contains x b

let rec depth = function
   | Leaf _ -> 0
   | Node (_, _, _, a, b) ->
       1 + max (depth a) (depth b)

The first is to test whether a key x is in the tree or not. The second is to calculate the depth of a tree.

Insert

For insertion, there is a known trick found by Chris Okasaki in this paper (Red-black trees in a functional setting). Okasaki is also the author of the famous book, Purely Functional Data Structures.

The idea of insertion is: always insert a Red node at some leaf node so that property 3 holds, then fix property 2 up to the root. There are four cases where property 2 could break. For all the four cases, the grand parent is always black, and there will be two red nodes right in the below. Please see the picture on page 3 of the paper.

Then the buttom-up fixing goes up the root. The root may be changed into Red. To hold property 1, simply change the color of the root to Black.

let insert x (t: RBTree<'a>) =
   let balance = function
       | Black, z, size, Node (Red, y, _, Node (Red, x, _, a, b), c), d
       | Black, z, size, Node (Red, x, _, a, Node (Red, y, _, b, c)), d
       | Black, x, size, a, Node (Red, z, _, Node (Red, y, _, b, c), d)
       | Black, x, size, a, Node (Red, y, _, b, Node (Red, z, _, c, d))->
           Node (Red, y, size, Node (Black, x, (getSize a b)+1, a, b), Node (Black, z, (getSize c d)+1, c, d))
       | color, b, _, c, d  -> // grandparent is red, does not change
           Node (color, b, (getSize c d)+1, c, d)

   let rec ins = function
       | Leaf _ ->
           Node (Red, x, 1, Leaf(Black), Leaf(Black))
       | Node (color, y, size, a, b) as s  ->
           if x < y then balance(color, y, (getSize a b)+2, ins a, b)
           elif x > y then balance(color, y, (getSize a b)+2, a, ins b)
           else s

   match ins t with
       | Node (Red, y, size, a, b)->
           Node (Black, y, size, a, b)
       | t -> t

Order statistics

Let’s have some rest before going to the monstrous delete operation. To get the nth smallest key in the tree, we can add a size field to each node, which is the size of the subtree rooted at the node. Thus the size of Node(_, _, _, left, right) is defined as

size = left.size + right.size + 1

and, as a boundary condition, a leaf node has size 0.

With this size field, the nth function is given below:

let nth n t =
   let rec nth' n = function
       | Leaf _->
           failwith "can't find the nth member"
       | Node (_, v, size, l, r)->
           let lsize=getSize l (Leaf(Black))
           if lsize+1=n then
               v
           elif lsize+1>n then
               nth' n l
           else
               nth' (n-lsize-1) r
   if n >= size t || n < 0 then
       failwith "nth out of range"
   else
       nth' (n+1) t

Delete

The idea of delete is simple (if we don’t care about the three properties of Red Black Trees): if the node x to be deleted has less than two children, then simply delete it, and lift its child (if any) to its position; if the node has two children, that means there must a node that is just bigger than x in its right child, replace x with that value, now the problem reduces to delete the smallest node in x’s right child.

First let’s write the function to get the smallest value in a tree:

let rec min = function
   | Leaf _->
       None
   | Node (_, value, _, Leaf(Black), _)->
       Some(value)
   | Node (_, _, _, l, r)->
       min l

Then we continue to code the delete logic:

        let rec del value = function
           | Leaf(color) ->
               failwith "RBT delete failed because x is not in the tree"
           | Node (color, y, _, a, b) as s ->
               if value < y then balance(color, y, 0, del value a, b)
               elif value > y then balance(color, y, 0, a, del value b)
               else
                   match (a, b) with
                   | (Leaf(Black), Leaf(Black)) -> // without non-leaf children
                       Leaf(Black++color)
                   | (Leaf(Black), Node (nc, nv, size, nl, nr))
                   | (Node (nc, nv, size, nl, nr), Leaf(Black)) ->  // with a single child
                       Node (color++nc, nv, size, nl, nr)
                   | (l, r) ->
                       //find the successor (smallest element in a right child),
                       //replace the key with the successor's
                       //and delete relevant node
                       match (min r) with
                       | Some(v) ->
                           balance(color, v, 0, l, del v r)
                       | None ->
                           failwith "impossible: can't find the successor"

The hard part is how to balance the tree after the deletion of the node. If we delete a black node, then Property 3 cannot hold anymore. So we must make rotations and do some rearrangements on colors to make Property 3 to hold again.

Think this problem in another way: after the deletion of a Black node and we connect its single child or leaf node into it, the node becomes darker. Conceptually Property 3 still holds if we think the node has two “Black”s in it. We mark this node as DoubleBlack and the task of balance is to eliminate this DoubleBlack node!

I defined an operator ++ to plus two colors:

let (++) color1 color2 =
   match color1,  color2 with
       | Red, Black
       | Black, Red ->
           Black
       | Black, Black ->
           DoubleBlack
       | _, _ ->
           Red

Let’s define the invariant for the balance function as follows: no node under pv is DoubleBlack.

Consider pv is the root of the subtree, and its left child is rooted at x, which is a DoubleBlack. Let’s consider all the cases for its right child. (After this is done, then we can work out when x is located as the right similarly.)

Case 1. If pv is black and its right child rv is red. Then do a left-rotation and make the DoubleBlack node x one level down (so that we reduce the problem):

(notations in the picture: nodes in circle are single nodes with colors, blue represents either red or black; a symbol represents a tree rooted at the symbol, for some cases the colors of the symbols are also noted.)

Case 2: If pv is of anycolor, and its right child is black. The three sub cases are shown below:

And the balance procedure is as follows:

let rec balance = function
   // invariant: the two children rooted at node pv have same "black" depth
   // but the root node itself may become a "double" black node
   | Black, pv, _, x, Node (Red, rv, size, ra, rb) when (isDoubleBlack x)-> // case 1
       balance(Black, rv, 0, balance(Red, pv, 0, x, ra), rb) // do left-rotate and continue balance
   | pc, pv, _, x, Node (Black, rv, size, ra, rb) when (isDoubleBlack x)-> // case 2.a 2.b and 2.c
       if isBlack ra && isBlack rb then // 2.a: reduce a black on both side
           let tempNode=Node (Red, rv, size, ra, rb)
           Node (pc++Black, pv, (getSize x tempNode)+1, blackify x, tempNode) // reduces the double node to root
       elif isBlack rb then // 2.b: do a right rotation in the right child, so rb becomes red and recudes to the "else" case
           match ra with
           | Node (_, rav, _, raa, rbb)->
               let tempNode1= Node (Red, rv, (getSize rbb rb)+1, rbb, rb)
               let tempNode2=Node (Black, rav, (getSize raa tempNode1)+1, raa, tempNode1)
               balance(pc, pv, 0, x, tempNode2)
           | _->
               failwith "impossible error"
       else // 3.c
           let tempNode=Node (Black, pv, (getSize x ra)+1, blackify x, ra)
           Node (pc, rv, (getSize tempNode rb)+1, tempNode, blackify rb)

   // when doubleblack x is on the right
   | Black, pv, _, Node (Red, lv, _, la, lb), x when (isDoubleBlack x)->
       balance(Black, lv, 0, la, balance(Red, pv, 0, lb, x))
   | pc, pv, _, Node (Black, lv, size, la, lb), x when (isDoubleBlack x)->
       if isBlack la && isBlack lb then
           let tempNode=Node (Red, lv, (getSize la lb)+1, la, lb)
           Node (pc++Black, pv, (getSize tempNode x)+1, tempNode, blackify x)
       elif isBlack la then
           match lb with
           | Node (_, _, lbv, lba, lbb)->
               let tempNode1=Node (Red, lv, (getSize la lb)+1, la, lba)
               let tempNode2=Node (Black, lbv, (getSize tempNode1 lbb)+1, tempNode1, lbb)
               balance(pc, pv, 0, tempNode2, x)
           | _->
               failwith "impossible error"
       else  
           let tempNode=Node (Black, pv, (getSize lb x)+1, lb, blackify x)
           Node (pc, lv, (getSize la tempNode)+1, blackify la, tempNode)   
       
   | a, b, _, c, d->
       Node (a, b, (getSize c d)+1, c, d)

Remember that this balance only balances the two children of any tree. Thus the root can still be DoubleDark (similar case in insert). We need a final check on the root:

match del x t with
| Node (_, y, size, a, b) ->
   Node (Black, y, size, a, b)
| Leaf _ -> Leaf(Black)

Validation for the implementation

It is important to have some test to verify the correctness of the above implementation. My idea is to do some random insertions and deletions, and then after each insertion or deletion, check whether the tree still holds the three properties.

The validate function is given below:

let validate t =
   let rec validate = function 
       | Leaf Black -> true, 0
       | Leaf DoubleBlack ->
           printfn "DoubleBlack in a Leaf node"
           false, 0
       | Leaf Red ->
           printfn "Red in a Leaf node"
           false, 0
       | Node (Black, _, _, a, b) ->
           let lr, ln = validate a
           let rr, rn = validate b
           lr && rr && (ln = rn), ln + 1
       | Node (Red, _, _, a, b) ->
           match a, b with
               | Node (Red, _, _, _, _), _
               | _, Node (Red, _, _, _, _) ->
                   let lr, ln = validate a
                   printfn "Red-Red in a tree"
                   false, ln
               | _ ->
                   let lr, ln = validate a
                   let rr, rn = validate b
                   lr && rr && (ln = rn), ln
       | Node (DoubleBlack, _, _, a, _) ->
           let lr, ln = validate a
           printfn "DoubleBlack in an internal node"
           false, ln
   match t with
   | Leaf Black -> true, 0
   | Leaf _ -> false, 0
   | Node (Black, _, _, _, _) ->
       validate t
   | Node (Red, _, _, _, _)
   | Node (DoubleBlack, _, _, _, _) ->
       printfn "root must be Black"
       false, 0

I did use the above test case to catch one bug in the implementation, the original code for the final deletion step is:

match del x t with

| Node (_, y, size, a, b) ->

Node (Black, y, size, a, b)

| t -> t

This code has a bug that when the tree is deleted to empty, its representation is Leaf (DoubleBlack) rather than Leaf(Black)!

I did a performance test against with F#’s set, which uses AVL tree:

    // Fisher-Yates shuffle
   let randperm n =
       let p = Array.init n (fun i -> i)
       let r = new System.Random(1)
       for i=n-1 downto 1 do
           let j = r.Next(i+1) // random number 0 <= j <= i
           let tmp = p.[j]
           p.[j] <- p.[i]
           p.[i] <- tmp
       p

   let lst = randperm 1000000 |> Seq.toList

   let rb = RBTree.ofList lst |> ignore

   open Microsoft.FSharp.Collections

   let av = Set.ofList lst |> ignore

   //timing on a AMD E-350 Netbook (1.6 GHz)
   //    Real: 00:00:44.635, CPU: 00:00:44.694, GC gen0: 218, gen1: 61, gen2: 4
   //
   //    val rb : unit = ()
   //
   //    >
   //    Real: 00:00:17.448, CPU: 00:00:18.766, GC gen0: 76, gen1: 20, gen2: 2
   //
   //    val av : unit = ()
   //

As we can see, for a large set with one million elements, my implementation costs about two times of that of AVL tree in F# Core library. That was not bad. I shall devote more time into performance and see if we can do more improvement. One optimization is to reduce redundant tree node allocation.

Tree traversal

Let’s talk some interesting stuff about tree traversal at last. I did implement some traversal functions:

let rec map f = function
   | Leaf _ as l -> l
   | Node (color, x, size, l, r) -> Node (color, f x, size, map f l, map f r)

let rec iter f = function
   | Leaf _ -> ()
   | Node (color, x, size, l, r) ->
       iter f l
       f x
       iter f r
  
let rec fold f t acc =
   match t with
   | Leaf _ -> acc
   | Node (color, x, size, l, r) -> fold f l (f x (fold f r acc))

They are straightforward recursive functions, and not tail-recursive. For a balanced tree like red black tree, tail-recursive is not important because the tree does not go too deep (a tree with 1 to 1000000 in it has a depth of 26). But it is just fun to make tem tail-recursive by using continuations. The standard accumulating variable trick does not work here because we have two recursive calls.

Let’s write the tail-recursive form of map function:

let map2 f t =
   let rec map f t k =
       match t with
       | Leaf _ as l -> k l
       | Node (color, x, size, l, r) ->
           map f l (fun lmap ->
               map f r (fun rmap -> k (Node(color, f x, size, lmap, rmap))))

   map f t (fun x -> x)

We can also use a Monad to abstract the continuation pattern out:

type Cont<'a,'r> = ('a -> 'r) -> 'r

type ContBuilder() =
   member x.Return (a):Cont<'a,'r> = fun k -> k a
   member x.Bind (c:Cont<'a,'r>, f:'a -> Cont<'b,_>) =
       (fun k -> c (fun a -> f a k))
   member this.Delay(f) = f()

let cont = ContBuilder()

let map3 f t =
   let rec map f = function
       | Leaf _ as t -> cont { return t }
       | Node (color, x, size, l, r) ->
           cont {
               let! lm = map f l
               let! rm = map f r
               return Node (color, f x, size, lm, rm)
           }
   map f t (fun x -> x)

A note on <'a when 'a:comparison>

In our definition of the Red Black Tree:

type RBTree<'a when 'a:comparison> =
   | Node of Color * 'a * int * RBTree<'a> * RBTree<'a>
   | Leaf of Color

The type signature has an ad hoc construct: when ‘a: comparison. Actually the only other constraint on ‘a is equality. F# only supports these two type constraints.

Type constraints behave like type classes in Haskell. :equality is Eq, :comparison is Ord. However, type classes in Haskell are much flexible.

Details could be found on this blog article by Don Syme.

Summary

In this article, I have implemented a functional Red Black Tree with dynamic order statistics. The tree is complete and has been well tested. Before any performance tuning, the running speed of it is within a small factor (2-3) of the AVL tree in F# Core library.

This article shows many important features of functional programming: 1) pattern matching and algebraic data types help implement complex data structures; 2) recursion is everywhere and the correctness of the program could be proved by inductive reasoning; and finally 3) a light introduction to tree traversal and continuation and Monad in F#.

N-queen problem, F#, Scheme and Haskell!

2011-08-28T20:30:00.001+08:00

I am busy with my research and other F# work, I haven’t written blog articles recently. But I always keep an eye on interesting things on F#. The F# version of “Write Yourself a Scheme in 48 hours” by Loca Bologonese really caught my eyes. Luca took the Haskell version and translated it into F# with the main structure remained the same. This gives me (an F# programmer) a great opportunity to learn Haskell by reading the F# and the Haskell programs side by side! It is really a good learning experience (previously I know a little Haskell).

N-Queen in F# and Scheme

The fun part was that I wished to write a functional version of n-queen program in Scheme and see the execution speeds of the F# evaluator and the Haskell evaluator!

n-queen problem is one of my favorite puzzles. I learned depth first search/backtracking with this problem. It is a good exercise for how to write a non-trivial recursive program other than factorial or Fibonacci sequence.

As the readers of this blog are probably F# users, I gave an F# version first in case you are not familiar with Scheme. The scheme version of the program uses exactly the same idea and same data structures.

The F# version:

// check if the two queens attack each other
let conflict (i1, j1) (i2, j2) = 
    if i1 = i2 || j1 = j2 || abs(i1-i2) = abs(j1-j2) then
        true
    else
        false
    
// check if a new position is valid to a list of other (valid) positions
let check (newPos:int*int) (posList: (int*int) list) = 
    posList
    |> List.exists (fun pos -> conflict newPos pos)
    |> not

// the backtracking procedure to get all solutions
let rec searchQueen i n sols = 
    if i = n then 
        sols
    else
        let newSols = 
            sols
            |> List.map (fun sol -> 
                [0..n-1]
                |> List.filter (fun j -> check (i,j) sol)
                |> List.map (fun j -> (i,j)::sol)
                )
            |> List.concat
        searchQueen (i+1) n newSols

let allSolutions = searchQueen 0 8 [[]]

allSolutions |> List.length

The scheme version:

;; helper functions
(define (accumulate op init seq)
(if (null? seq)
      init
      (op (car seq)
          (accumulate op init (cdr seq)))))

(define (accumulate op init seq) (if (null? seq) init (op (car seq) (accumulate op init (cdr seq)))))

(define (flatmap proc seq)
(accumulate append '() (map proc seq)))

(define (enumerate-interval low high)
(if (> low high)
'()
(cons low (enumerate-interval (+ low 1) high))))

;; data structures and basic functions for board
(define (make-queen row col) (list row col))
(define (get-row queen) (car queen))
(define (get-col queen) (car (cdr queen)))

(define (same-col? q1 q2) (= (get-col q1) (get-col q2)))
(define (same-diag? q1 q2)
(=
(abs (- (get-row q1) (get-row q2)))
(abs (- (get-col q1) (get-col q2)))))

(define (attacks? q1 q2)
(or (same-col? q1 q2) (same-diag? q1 q2)))

(define (safe? newq qlist)
(cond
    ((null? qlist) #t)
    ((attacks? newq (car qlist)) #f)
    (else (safe? newq (cdr qlist)))))

(define (safe-board? qlist)
(let
    ((newq (car qlist))
     (rest (cdr qlist)))
    (safe? newq rest)))

;; the depth-first search for queens
(define (queens board-size)
(define (queen-rows k sols)
    (if (= k board-size)
        sols
    (queen-rows (+ k 1)
          (filter
            (lambda (board) (safe-board? board))
            (flatmap
          (lambda (rest-queens)
            (map (lambda (new-col)
               (cons (list k new-col) rest-queens))
             (enumerate-interval 1 board-size)))
          sols)))))
(queen-rows 0 (list '())))

(length (queens 8))

Here are some correspondences between the F# version and the scheme version:

1. In F#, each queen position is represented as a tuple (int * int), while in Scheme, each queen position is a list with length 2.

2. In the search procedure, the main logic is: given a list of sub-solutions (sols), expand every one with all possible columns for the new position. This gives a list of lists of lists of positions. There should be a flatten function to transform it into a list of lists of positions. In F#, it is List.concat, in Scheme, it is flatmap.

Run n-queen in the Scheme interpreters

The fun part is to run the n-queen program in the F# version interpreter. However, throwing the above scheme program into the interpreter gives errors. (It runs well in CzScheme)

Here are the limitations that I encountered:

1. The interpreter only supports one-line definitions. Overcoming this limitation is easy – just write the functions in one-lines.

2. The interpreter does not support (let). let syntax gives a word for some expression. The easy way is to use the expression directly.

3. Some built-in functions are not included in the interpreter, e.g. abs and append. Implementing them is quite easy.

4. The interpreter uses Console.ReadLine() to read from the input, which only supports 256 characters per line by default. The following lines reset the buffer size to allow more characters per line:

        let inputBuffer = Array.create 1024 0y          let inputStream = Console.OpenStandardInput(inputBuffer.Length)         Console.SetIn(new IO.StreamReader(inputStream))

5. Does not support nested-define. Nested-define could be implemented as let + lambda. But I didn’t solve this problem in this post. I just remove nested-defines.

Here is the runnable Scheme program:

(define (abs x) (if (negative? x) (- 0 x) x))
(define (append list1 list2) (if (null? list1) list2 (cons (car list1) (append (cdr list1) list2))))
(define (accumulate op init seq) (if (null? seq) init (op (car seq) (accumulate op init (cdr seq)))))
(define (flatmap proc seq) (accumulate append (quote ()) (map proc seq)))
(define (enumerate-interval low high) (if (> low high)    '()      (cons low (enumerate-interval (+ low 1) high))))
(define (make-queen row col) (list row col))
(define (get-row queen) (car queen))
(define (get-col queen) (car (cdr queen)))
(define (same-col? q1 q2) (= (get-col q1) (get-col q2)))
(define (same-diag? q1 q2)   (=    (abs (- (get-row q1) (get-row q2)))    (abs (- (get-col q1) (get-col q2)))))
(define (attacks? q1 q2) (or (same-col? q1 q2) (same-diag? q1 q2)))
(define (safe? newq qlist) (if     (null? qlist) #t    (if (attacks? newq (car qlist)) #f (safe? newq (cdr qlist)))))
(define (safe-board? qlist) (safe? (car qlist) (cdr qlist)))
(define board-size 7)
(define (queen-rows k sols) (if (= k board-size)   sols (queen-rows (+ k 1)           (filter             (lambda (board) (safe-board? board))            (flatmap        (lambda (rest-queens)        (map (lambda (new-col)      (cons (list k new-col) rest-queens))       (enumerate-interval 1 board-size)))       sols)))))
(length (queen-rows 0 (list (quote ()))))

The above program runs OK in both F# and Haskell interpreter at about the same speed. However, when board-size equals 8, F# version gets a Stackoverflow runtime error.

I think it is enough for a post, I will see if I can solve these problems

Extracting top words from titles and abstractions in MIX11 presentations

2011-04-18T11:32:00.001+08:00

The videos of MIX11 (Apr. 11-14) conference are available at http://channel9.msdn.com/Events/MIX/MIX11. To get a sense of what’s going on recently, I wrote an F# script to count the top occurring words in the titles and abstractions of all the presentations.

Title statistics:

"windows" : 41
"phone" : 26
"azure" : 14
"web" : 14
"html5" : 11
"silverlight" : 9
"net" : 9
"7" : 8
"applications" : 8
"data" : 8
"platform" : 7
"application" : 7
"javascript" : 7
"new" : 7
"ux" : 6
"what’s" : 6
"asp" : 6
"boot" : 5
"camp" : 5
"building" : 5

From the statistics of the words in titles, we can find that windows mobile phone and azure cloud platform are the hottest topics. HTML5, the next standard of web page technology, has been always been a focus by Microsoft. Silverlight still has its heat. It is a good company to HTML5, I think it will have many applications in in-house web applications; while HTML5 has more support and available in all browsers across platforms.

Abstract statistics:

"web" : 62
"session" : 55
"windows" : 49
"new" : 40
"applications" : 36
"phone" : 36
"learn" : 32
"we’ll" : 29
"use" : 25
"using" : 24
"silverlight" : 23
"data" : 23
"net" : 22
"javascript" : 21
"azure" : 21
"developers" : 20
"come" : 17
"microsoft" : 17
"features" : 17
"one" : 16

and the whole F# program is actually short -- only 40 lines! You have everything there: download webpages, get the titles and titles, stopword removing and word counting, sorting...

Code Snippet
 open System
open System.Net
open System.Text.RegularExpressions
 
let fetchUrlSimple (url:string) = 
    let req = WebRequest.Create(url) 
    let response = req.GetResponse()
    use stream = response.GetResponseStream()
    use streamreader = new System.IO.StreamReader(stream)
    streamreader.ReadToEnd()
 
let topKWords (docs:string seq) K = 
    let separator = [|' '; '\r'; '\n'; '-'; '.'; ',' ; '\t'; '!'; '?'; '\''; ';'; '/' |]
    let stopwords = Set("a,able,about,across,after,all,almost,also,am,among,an,and,any,are,as,at,be,because,been,but,by,can,cannot,could,dear,did,do,does,either,else,ever,every,for,from,get,got,had,has,have,he,her,hers,him,his,how,however,i,if,in,into,is,it,its,just,least,let,like,likely,may,me,might,most,must,my,neither,no,nor,not,of,off,often,on,only,or,other,our,own,rather,said,say,says,she,should,since,so,some,than,that,the,their,them,then,there,these,they,this,tis,to,too,twas,us,wants,was,we,were,what,when,where,which,while,who,whom,why,will,with,would,yet,you,your".Split ',')
 
    docs
    |> Seq.map (fun doc ->
        doc.Split(separator, StringSplitOptions.RemoveEmptyEntries)
        |> Seq.map (fun word -> word.ToLower())
        |> Seq.filter (fun word -> not (stopwords.Contains(word)))
        )
    |> Seq.concat
    |> Seq.groupBy (fun x->x)
    |> Seq.map (fun (word, wordSeq) -> (word, wordSeq |> Seq.length))
    |> Seq.sortBy (fun (_, wordCnt) -> - wordCnt)
    |> Seq.take K
    |> Seq.toList
 
let rawpage = fetchUrlSimple @"http://channel9.msdn.com/Events/MIX/MIX11" 
let titles, abstracts =
    let page = Regex.Replace(rawpage, "&#?[a-z0-9]+;", " ")
    let titleMatches = Regex.Matches(page, "class=\"title\">(.*?)</a>")
    let abstractMatches = Regex.Matches(page.Replace("\n"," "), "class=\"description\">(.*?)</div>")
    let matchesToSeq (matches: MatchCollection) = 
        seq {
            for m in matches do
                yield m.Groups.[1].Value
        }
    matchesToSeq titleMatches, matchesToSeq abstractMatches
 
 
topKWords titles 20
topKWords abstracts 20
 

My colleague Defu Lian wrote a C# version for topKWords function using LINQ:

Code Snippet
 static List<Tuple<string, int>> topKWords(IEnumerable<string> docs,int K)
{
    char[] separator = { ' ', '\r', '\n', '-', '.', ',', '\t', '!', '?', '\'', ';', '/' };
    var stopwords = new HashSet<string>("a,able,about,across,after,all,almost,also,am,among,an,and,any,are,as,at,be,because,been,but,by,can,cannot,could,dear,did,do,does,either,else,ever,every,for,from,get,got,had,has,have,he,her,hers,him,his,how,however,i,if,in,into,is,it,its,just,least,let,like,likely,may,me,might,most,must,my,neither,no,nor,not,of,off,often,on,only,or,other,our,own,rather,said,say,says,she,should,since,so,some,than,that,the,their,them,then,there,these,they,this,tis,to,too,twas,us,wants,was,we,were,what,when,where,which,while,who,whom,why,will,with,would,yet,you,your".Split(','));
    return docs.SelectMany(doc => doc.Split(separator, StringSplitOptions.RemoveEmptyEntries)
                        .Select(word => word.ToLower())
                        .Where(word => !stopwords.Contains(word)))
        .GroupBy(a => a)
        .Select(wordgroup => Tuple.Create<string,int>(wordgroup.Key, wordgroup.Count()))
        .OrderByDescending(w2c => w2c.Item2)
        .Take(K)
        .ToList();
}
 

A Note on F# Quotations

2011-04-10T20:44:00.001+08:00

I worked on extending the F# ODSL (Optimization Domain Specific Language, link1 and link2) during the weekend. This is the first time I use F# quotations in a non-trivial fashion. I tried some quotation examples in Programming F# when I was learning the language and read code examples in F# PowerPack (for LINQ integration) and other libraries before.

Only after I actually write some programs in it, I can appreciate this F# feature more and have some of my own thoughts. In this post, I’d like to share these thoughts. I don’t intend to go into detailed F# code, but maybe in the future I will write a separate blog on my work on extending the F# ODSL.

General idea of F# Quotations

The great idea of quotation at least traces back to Lisp, where program is also a kind of data – the execution behavior of a piece of program is completely controllable by the user, just treat it as input data and write a custom evaluator for it. The default Lisp evaluator is eval, we can easily write a custom one to change the default behavior [LispEval].

In F#, we can also treat a piece of F# code as data by quoting it:

let quotedProgramAsData = 
    <@
        let y = 1 + 2 * 30
        20 * y
        // code here:
        // an F# program using a subset of F# language features
    @>

When we create the variable quotedProgramAsData, the .Net runtime will not compute the expression in the quotation; instead it generates the following value:

val quotedProgramAsData : Quotations.Expr<int> =

val quotedProgramAsData : Quotations.Expr<int> =
  Let (y,
     Call (None, Int32 op_Addition[Int32,Int32,Int32](Int32, Int32),
           [Value (1),
            Call (None, Int32 op_Multiply[Int32,Int32,Int32](Int32, Int32),
                  [Value (2), Value (30)])]),
     Call (None, Int32 op_Multiply[Int32,Int32,Int32](Int32, Int32),
           [Value (20), y]))

Let’s visualize it in a tree form:

F# compiler generates this tree structure by free and the rest is how you deal/evaluate this tree. If the above code is not in quotation, F# compiler will generate code that put the value of the second parameter (1 + 2 * 30) of Let to the first parameter (y), and continue to generate code for the third parameter (20 * y). But because it is in the quotation, only the expression tree is generated without any explicit execution behavior for them.

With different purposes, we can write different evaluators for the quoted F# code. In the following I list some of the application areas.

Domain specific language (DSL)

There are different kinds of domain specific languages. In his book Domain Specific Languages, Martin Fowler summarizes them into external and internal ones:

DSLs come in two main forms: external and internal. An external DSL is a language that's parsed independently of the host general purpose language: good examples include regular expressions and CSS. External DSLs have a strong tradition in the Unix community. Internal DSLs are a particular form of API in a host general purpose language, often referred to as a fluent interface. The way mocking libraries, such as JMock, define expectations for tests are good examples of this, as are many of the mechanisms used by Ruby on Rails. Internal DSLs also have a long tradition of usage, particularly in the Lisp community. (http://martinfowler.com/books.html#dsl)

Quotations can be used as a very good language feature for implementing internal domain specific languages.

As stated as before, the quoted program is anyway an F# program. How could it be look like a domain specific one? Remember F# has a rich set of syntax while a domain language takes a small subset of it is usually enough expressive.

Take a look at the following program, which expresses the quadratic optimization for Support Vector Machines:

// solve the SVM using ODSL
let dsl_solver = 
    let index = [ 0..n-1]
    Solver 
        <@
            qp()
            let alpha = vararray1 (index)
            maximise (sum index (fun i -> alpha.[i]) - 
                        (sum index (fun i -> sum index (fun j -> 
                            coef.[i,j] * alpha.[i] * alpha.[j])))) // Eq. (1)
            where 
                [
                    foreach index (fun i -> 0.0 <= alpha.[i]) // Eq. (2)
                    foreach index (fun i -> alpha.[i] <= C) // Eq. (2)
                    sum index (fun i -> alpha.[i] * y.[i]) = 0. // Eq. (3)
                ]
        @>

and its optimization formulation:

Subject to constraints:

The correspondence between the DSL code and mathematical formulates is very clear. In the future work, we can even vectorize this DSL to eliminate the usage of sum: range -> lambda (int->value) -> value function, which will make the program and the formulas more similar.

Note that in the above example qp, maximise, sum and where are not magic, but are all F# functions whose behaviors are defined in the DSL. For example, F# function qp() only indicates the following code is a quadratic programming.

Compare this code with the code in my SVM post: it is clearer and thus easier to write.

The restriction of such a DSL inside F# is that the syntax of this DSL should follow that of F# -- it should be a valid program recognizable by the F# lexer and parser. This is why DSL people like Lisp which has so simple and flexible syntax (just brackets), but sometimes it is also boring and confusing when one gets lost in the brackets.

High performance computing

We know that F#’s Seq module and C#’s LINQ share some common features, e.g. chained operations and the laziness. Because these operations are so fundamental, LINQ team has spent enormous time optimizing LINQ; on the other hand, F#’s implementation is quite a standard one without heavy optimization because F# compiler team is not big and they have other important tasks to do. (F#’s Seq module actually does some optimization, e.g. for arrays and lists special routines are called instead of the general ones for IEnumerable<T> objects. But LINQ does more!)

So one idea of speeding up sequence operations in F# is to use LINQ’s equivalent functions. This could be implemented by putting F#’s sequence operations in a quotation and rewrite them in LINQ expression as done in F# PowerPack. Interested readers can read the source code of LINQ for F# in PowerPack.

The above example talks about basic data structures. As a data mining guy myself, let’s move to numerical computing: we can quote a piece of numerical code in F# and translates it into Fortran/C and a JIT Fortran/C compiler compiles the code translated code into native machine code which uses special instruction set in that platform (e.g. vectorized instructions in P4 CPU). Or even we can compile this piece of F# into GPU and utilize the parallel computing there. [Syme_ML06]

But any fancy optimization has some overhead too. For example, the compiling time may cost some time.

Education: Compiler courses

The last application area is education. There are few CS departments teaching functional programming. Some of them will mention it when teaching programming language concepts. F# is a push for FP into industry and main stream. Hope in the near future, schools will open tiny or selected courses for F#. See the recent F# in Education Workshop [Education] for details.

Ok. Let’s focus on compiler courses. Compilers, different from functional programming, are taught in nearly all CS departments. However the projects in compiler courses are a huge pain to students. They are so complicated! And the focus of the first half of the course (lexing and parsing) and the second half (code generation and optimization) are kind of separate. The theories behind the two halves are different. Sometimes, the course instructor focuses too more on lexing and parsing and don’t have enough time for students to work on code generation and optimization, which in my mind are more important in a compiler course. While lexing and parsing will occur in other CS courses too, e.g. computational theory, code generation and optimization will only be in a compiler course.

By using the F# quotations, F# compiler generates the expression tree for free. Just so convenient and straightforward! Considering that the F# syntax you can put in quotation is a capable set of imperative language + functional language, the expression tree for us to do the code generation and optimization is non-trivial. By using tools like .Net Reflector, we can also study how the standard F# compiler generates the same piece of code and learn the tricks there. Students can work on code generation for closures, and tail-recursive-calls – important language features that are not implemented or fully implemented in many main stream languages.

References

Some references in the above text:

[LispEval] Christian Queinnec , Lisp in Small Pieces.

[Syme_ML06] Leveraging .NET Meta-programming Components from F#, ML Workshop, 2006.

[Education] F# in Education Workshop, http://research.microsoft.com/en-us/events/fsharpined/.

A very good paper by Leijen and Meijer:

[Leijen&Meijer] Domain specific embedded compilers, Sigplan Notices, vol. 35, no. 1, pp. 109-122, 2000.

F# quotation extensively uses active patterns for expression tree pattern matching:

[Syme_Active] Extensible Pattern Matching via a Lightweight Language Extension, ICFP 2007.

Here are some relevant blog posts:

[Petricek] F# Overview (IV.) - Language Oriented Programming, http://tomasp.net/blog/fsharp-iv-lang.aspx.

[RubyDSL] Building a DSL in Ruby, http://jroller.com/rolsen/entry/building_a_dsl_in_ruby.

Two nice answers by Petricek & Harrop discussing Quotations (F#, Ocaml and Lisp) and DSL on Stackoverflow:

http://stackoverflow.com/questions/4317912/confusing-about-f-quotations-and-pattern-matching-in-meta-programming.

Acknowledgement

I’d like to thank Nathan Brixius for discussions on F# ODSL!

Support vector machines (SVMs) in F# using Microsoft Solver Foundation

2011-02-24T15:40:00.001+08:00

Support vector machines are a super star in machine learning and data mining in the past decade. It reheats statistical learning in machine learning community. It is also one of the classifiers that work practically with applications to other disciplines such as bioinformatics.

In this post, I will give a brief review of the theory of this classifier and show an F# implementation using the quadratic programming (QP) solver in Microsoft Solver Foundation (MSF).

It may seem dull to start with mathematics. But the math in SVMs is actually quite intuitive and would be quite easy if one has some training in linear algebra and calculus. Also understanding the intuition behind the formulae is more important in the formulation itself. There are quite a lot of good SVM tutorials written by researchers. Thus I only write the important formulas that occur in the implementation. The reader could link the F# code with formulas directly.

My note is based on the machine learning course by Prof. Jaakkola [Jaakkola] and the Matlab implementation therein. The course material is highly comprehensive yet accurate.

Background review for SVMs

Given a binary labeled dataset, where and is 1 or -1. Think about the data points plotted in a d dimensional (Euclidean) space, the linear SVM classifier is a hyperplane in the space and it “best” separates the two kinds of data points. Here “best” means that the hyperplane should have the largest margin, which is the distance from the plane to the sides of labeled points. Intuitively the larger the margin is, the more robust and confident the classifier is. As if the hyperplane shakes a little, it still classifies well for its being far from both sides.

Let’s take d=2 as the discussing example: data points are plotted in a 2-D plane. The aim is to draw a line to separate the two kinds of points such that the margin of the separation is maximized. In the following Figure, the line with maximum margin is shown. After some geometry work, the margin is calculated as , thus minimizing with the constraints that the two kinds of points are well separated () gives the max-margin hyperplane.

(Wikipedia picture)

Two important extensions to the above simplest SVM are:

1. Allowing imperfect separation, that is when a hyperplane could not separate the points perfectly, it is allowed to mis-separate some data points. This introduces a slack variable for each data point , when the point is on the correct side of the plane, the slack variable is 0, otherwise it is measures the distance it goes from the plane minus 1, if the values goes to below zero set it zero. (Please read the Note and reference section for more explanation.)

2. Kernel. The simplest SVM is a linear classifier in the original space. However, there are situations where linear classifier is not sufficient (even the best), one strategy is to do a non-linear transformation or mapping that maps a data point to another space, so that a linear classifier in that space, which is non-linear in the original space, does a better separation of the data points. The kernel trick is that we don’t need to define explicitly; only the definition of inner product of that space is required: .

With the two extensions, the new maximum margin objective becomes:

subject to

The dual form:

subject to

The data points with its alpha value greater than 0 are called support vectors.

With the in the dual problem solved, the SVM classification hyperplane is recovered by:

The threshold parameter could be calculated by substituting back to the support vectors:

Any and (So that make sure data point is on the margin boundary, not in-between) would give a , to stabilize the solution, it is common practice to take the median of off the possible .

The implementation

The solver is actually an Interior Point Solver, which could solve linear and quadratic programming with linear constraints.

Page 14-17 of MSF-SolverProgrammingPrimer.pdf show an example usage of the QP solver in Microsoft Solver Foundation. But the example is not complete and it does not demonstrate one very subtle part: how to set the coefficients for the quadratic terms.

The following code example shows a simple example with the the comment regarding the coefficient setting:

Example usage of the QP solver
 #r @"C:\Program Files (x86)\Sho 2.0 for .NET 4\packages\Optimizer\Microsoft.Solver.Foundation.dll"
open Microsoft.SolverFoundation.Common
open Microsoft.SolverFoundation.Solvers
open Microsoft.SolverFoundation.Services
 
 
// test QP
module TEST = 
(* minimize x^2 + y^2 + 3xy + 2x + y *)
(* notice that in the InteriorPointSolver,
   the coefficients for xy & yx should be the same (so only set ONCE!)
   if we set both 3xy and 0yx, the solver takes the later coef. 
*)
    let solver = new InteriorPointSolver()
    let _, goal = solver.AddRow("dual objective value")
    solver.AddGoal(goal, 0, true)
 
    let _, x = solver.AddVariable("x")
    let _, y = solver.AddVariable("y")
 
    solver.SetCoefficient(goal, x, Rational.op_Implicit(2))
    solver.SetCoefficient(goal, y, Rational.op_Implicit(1))
    // for terms like x-y (where x != y), set its coef for only one time!
    solver.SetCoefficient(goal, Rational.op_Implicit(3), x, y)
    //solver.SetCoefficient(goal, Rational.Zero, y, x)
    solver.SetCoefficient(goal, Rational.op_Implicit(1), x, x)
    solver.SetCoefficient(goal, Rational.op_Implicit(1), y, y)
    
    
    let param = new InteriorPointSolverParams()
    solver.Solve(param) |> ignore
 
    //solver.Result
    let objectiveValue = solver.GetValue(0).ToDouble()
 
    let x0 = solver.GetValue(1).ToDouble()
    let y0 = solver.GetValue(2).ToDouble()
    x0*x0 + y0*y0 + 3.0 * x0 * y0 + 2. * x0  + y0
    //x0*x0 + y0*y0 + 0.0 * x0 * y0 + 2. * x0  + y0
 

The implementation of SVM is a straightforward translation of equations (1) to (5). The following shows the svm learning (bulidSvm) and testing (svmClassify) functions:

SVM Implementation
 open System
 
 
type dataset = 
    { features: float array array; // (instance = float array) array
      mutable labels: int array; // 
    }
    with
    member x.NumSamples = x.features.Length
 
 
module Array = 
    let median (a:'a array) = 
        let sorted = Array.sort a
        sorted.[sorted.Length / 2]
        
module Kernel = 
    let linear a b = 
        Array.fold2 (fun acc p q -> acc + p * q) 0.0 a b
 
    let polynomial k a b = 
        let dot = linear a b
        Math.Pow(1.0 + dot, k |> float)
 
    let gaussian beta a b = 
        let diff = Array.fold2 (fun acc p q -> acc + (p-q)*(p-q)) 0.0 a b
        exp (-0.5 * beta * diff)
 
module SVM = 
    type svmmodel = {
        SVs:dataset;
        alpha:float array;
        kernel: float[] -> float[] -> float;
        w0:float;
        }
        with
        member x.NumSupporVectors = x.SVs.features.Length
 
 
    let buildSVM (ds:dataset) (C:float) (kernel:float[] -> float[] -> float) = 
 
        let n = ds.features.Length
        let C = Rational.op_Implicit(C)
        let zero = Rational.Zero
 
        // create a interior point solver, which solves the QP problem
        let solver = new InteriorPointSolver()
 
        // set the objective value / goal
        let _, goal = solver.AddRow("dual objective value")
 
        // false == maximizing the objective value
        // the value of goal is (1)
        solver.AddGoal(goal, 0, false) |> ignore
 
        // add the Lagangian variables \alpha_i and set their bounds (0 <= \alpha_i <= C)
        let alpha = Array.create n 0
        for i=0 to n-1 do 
            let _, out = solver.AddVariable("alpha_"+i.ToString())
            alpha.[i] <- out
            solver.SetBounds(out, zero, C)
 
        // add contraint: \sum_i \alpha_i * y_i = 0
        // equation (2)
        let _, sumConstraint = solver.AddRow("SumConstraint")
        solver.SetBounds(sumConstraint, zero, zero);
 
        for i=0 to n-1 do
            // set the coefs for the sum constraint
            // equation (2)
            solver.SetCoefficient(sumConstraint, alpha.[i], Rational.op_Implicit(ds.labels.[i]))
 
            // add the \alpha_i terms into the objective
            solver.SetCoefficient(goal, alpha.[i], Rational.One)
 
            // add the qudratic terms
            for j=0 to i do 
                // coef = y_i * y_j * K(x_i, x_j)
                let coef = float(ds.labels.[i] * ds.labels.[j]) * (kernel ds.features.[i] ds.features.[j])
                if i=j then
                    solver.SetCoefficient(goal, Rational.op_Implicit(-0.5 * coef), alpha.[i], alpha.[j])
                else
                    solver.SetCoefficient(goal, Rational.op_Implicit(-coef), alpha.[i], alpha.[j])
                
        // use the default parameters 
        let param = new InteriorPointSolverParams()
 
        
        solver.Solve(param) |> ignore
 
        // get the alpha values out
        let alphaValue = Array.init n (fun i -> solver.GetValue(i+1))
 
    (* print optimization result
        printfn "goal value = %A" (solver.GetValue(0).ToDouble())
        for i=1 to n do
            printfn "%A" (solver.GetValue(i).ToDouble())
    *)
        let alphaNew = new ResizeArray<Rational>()
        // extract the non-zero alpha values out and their corresponding support vectors
 
        
        let SVs = 
            let feats = new ResizeArray<float[]>()
            let labs = new ResizeArray<int>()
            let maxAlpha = Array.max alphaValue
            let threshold = maxAlpha * Rational.op_Implicit(1e-8)
            for i=0 to n-1 do 
                if alphaValue.[i] > threshold then
                    feats.Add(ds.features.[i])
                    labs.Add(ds.labels.[i])
                    alphaNew.Add(alphaValue.[i])
 
            { features = feats |> Seq.toArray;
              labels = labs |> Seq.toArray;
            }
 
        // solve w_0 in the primal form
    
        let alphaNZ = alphaNew |> Seq.toArray
 
        // equation (5)
        let w0 = 
            alphaNZ
            |> Array.mapi (fun i a -> 
                if a = C then 
                    None
                else 
                    let mutable tmp = 0.0
                    for j=0 to SVs.NumSamples-1 do
                        tmp <- tmp + alphaNZ.[j].ToDouble() * (SVs.labels.[j] |> float) * (kernel SVs.features.[i] SVs.features.[j])
                    Some ((float SVs.labels.[i]) - tmp)
                )
            |> Array.filter (fun v -> match v with None -> false | _ -> true)
            |> Array.map (fun v -> match v with Some v -> v | _ -> 0.0)
            |> Array.median
        
 
        // construct an svm record
        {
            SVs = SVs;
            alpha = alphaNZ |> Array.map (fun v -> v.ToDouble());
            kernel = kernel;
            w0 = w0;
        }
        
 
    let svmClassify (model:svmmodel) (ds:dataset) = 
        // equation (4)
        let vals = 
            ds.features
            |> Array.map (fun x ->
                let mutable sum = 0.0
                for i=0 to model.NumSupporVectors-1 do
                    sum <- sum + model.alpha.[i] * (float model.SVs.labels.[i]) * (model.kernel model.SVs.features.[i] x)
                sum + model.w0
                )
        let nCorrect = 
            Array.map2 (fun value label -> if (value > 0.0) && (label = 1) || (value < 0.0) && (label = -1) then 1 else 0) vals ds.labels
            |> Array.sum
    
        (float nCorrect) / (float ds.NumSamples), vals
 

Try on the iris data set we used in the Logistic Regression post:

let svm = buildSVM iris 10.0 Kernel.linear
let classifyResult = svmClassify svm iris

val svm : svm =
{SVs =
{features =
[|[|1.0; 6.2; 2.2; 4.5; 1.5|]; [|1.0; 5.9; 3.2; 4.8; 1.8|];
[|1.0; 6.3; 2.5; 4.9; 1.5|]; [|1.0; 6.7; 3.0; 5.0; 1.7|];
[|1.0; 6.0; 2.7; 5.1; 1.6|]; [|1.0; 5.4; 3.0; 4.5; 1.5|];
[|1.0; 4.9; 2.5; 4.5; 1.7|]; [|1.0; 6.0; 2.2; 5.0; 1.5|];
[|1.0; 6.3; 2.7; 4.9; 1.8|]; [|1.0; 6.2; 2.8; 4.8; 1.8|];
[|1.0; 6.1; 3.0; 4.9; 1.8|]; [|1.0; 6.3; 2.8; 5.1; 1.5|];
[|1.0; 6.0; 3.0; 4.8; 1.8|]|];
labels = [|-1; -1; -1; -1; -1; -1; 1; 1; 1; 1; 1; 1; 1|];};
alpha =
[|6.72796421; 10.0; 10.0; 10.0; 10.0; 6.475497298; 1.490719712;
8.547262902; 3.165478894; 10.0; 10.0; 10.0; 10.0|];
kernel = <fun:svm@226-32>;
w0 = -13.63716815;}

>
Real: 00:00:00.002, CPU: 00:00:00.000, GC gen0: 0, gen1: 0, gen2: 0

val classifyResult : float * float [] =
(0.97,
[|-2.787610618; -2.380530973; -1.42477876; -2.929203541; -1.681415929;
-1.96460177; -1.24778761; -6.106194692; -2.761061946; -2.973451328;
…. 4.203539824; 1.82300885|])

Notes and references

The optimization in SVMs could be treated as a general QP problem as shown in our implementation. However, when the number of the data points is big, the QP grows too big to handle by the QP solver.

As a special QP problem, SVM has a lot of novel solutions proposed in the machine learning research area, e.g. the SMO approach [Platt] in LibSVM implementation and the Ball Vector Machine approach [Tsang] transforming the special QP into a (discrete) computational geometry problem, minimum enclosing ball. For practical usage of SVM, one usually uses these dedicated approaches.

Besides the maximum margin interpretation, SVMs also have theory roots in functional regularization theory, in which the optimization

is viewed as minimizing the error term while controlling the complexity of the function and C is a tradeoff parameter between these two. The reader is referred to the teaching notes and papers by Prof. Poggio [Poggio].

[Burges] A Tutorial on Support Vector Machines for Pattern Recognition. http://research.microsoft.com/en-us/um/people/cburges/papers/svmtutorial.pdf

[Jaakkola] http://courses.csail.mit.edu/6.867/.

[Platt] Sequential Minimal Optimization: A Fast Algorithm for Training Support Vector Machines. ftp://ftp.research.microsoft.com/pub/tr/tr-98-14.pdf

[Poggio] http://www.mit.edu/~9.520/spring10/.

[Tsang] http://www.cs.ust.hk/~ivor/cvm.html.

Logistic Regression in F# using Microsoft Solver Foundation

2011-02-11T16:32:00.001+08:00

Logistic regression is a workhorse in data mining. Like decision trees and SVMs, it is a very standard classifier. If you have a labeled data, logistic regression definitely is one of the classifiers that should tried.

In this post, I’d like to show how to implement a logistic regression using Microsoft Solver Foundation in F#. Because we use the optimization in Microsoft Solver Foundation, the core part of a logistic regression contains only 20 – 30 lines of code.

If you are only interested in F# and L-BFGS, but not logistic regression, please only read the second section: L-BFGS solver in Microsoft Solver Foundation.

Background for logistic regression

Logistic regression model predicts the probability that a data instance x being labeled as y

A data instance is a vector with each dimension representing a value for some feature. For example, the Iris flower data set we are going to use in this post has 4 features: Sepal Length, Sepal Width, Petal Length and Petal Width. Each instance has a label y indicating the category. The above model only models binary labels. Logistic regression also has multiple categories version, see [Mitchell]. In this post, we focus on the binary version; however multiple versions should be straightforward once we know how to implement the binary one.

The model contains two parameters weight vector w and offset b, which should be estimated/learned on the training data set. In our derivations below, we will drop the offset b, because we can always add one feature (constant 1) to feature vector x’=[x 1] and hide b into the extended w’=[w b]. Sometimes, a model without the offset parameter could also perform equally well. Because equation (1) measures the probability, we could solve the model parameters by optimizing the joint probability estimated on dataset X.

If we assume that all the instances are independent and identically distributed (i.i.d.). Then the regularized joint probability is

where is the regularization term for avoiding over-fitting on the training data.

The log version of this likelihood is

We want to find the model parameter w that maximizes this likelihood. P(X) is a concave function with respect to w, so it has a global maximal.

What left is how to optimize this function. BFGS is a Quasi-Newton method that only needs first order gradient while a Newton family method needs second order explicitly. Microsoft Solver Foundation has a limited memory version of BFGS (L-BFGS). (L-BFGS is more popular than BFGS, sometimes people mention BFGS to mean L-BFGS.) There is a memory parameter m in L-BFGS, the space complexity of this algorithm is O(md), where d is the dimension of the variable space. Say, if you want to optimize a function that has 1,000,000 free variables, and set m to 10, then you need only a little more than 1M*10*8/1024^2=76M bytes.

The input for an L-BFGS optimizer is: 1) the function evaluator L(w) and 2) the gradient evaluator :

The output is the best that gives the maximum function value.

L-BFGS solver in Microsoft Solver Foundation

Microsoft Solver Foundation already implements an L-BFGS solver. The reader is suggested to read pages 24-27 & 72-75 of MSF-SolverProgrammingPrimer.pdf (on my machine: C:\Program Files\Microsoft Solver Foundation\3.0.1.10599\Documents) for detailed reference. Because the code for optimizing Rosenbrock function is in C#, here I translated it to F#. Note how I translate between delegate type System.Func<> and F# functions.

Btw, I find that using the Microsoft.Solver.Foundatoin.dll in Sho would be more convenient as it packages everything into one DLL.

#r @"C:\Program Files (x86)\Sho 2.0 for .NET 4\packages\Optimizer\Microsoft.Solver.Foundation.dll"
openMicrosoft.SolverFoundation.Common
openMicrosoft.SolverFoundation.Solvers
openMicrosoft.SolverFoundation.Services

letsolverParams = newCompactQuasiNewtonSolverParams()
letsolver = newCompactQuasiNewtonSolver()
//add variables
let_, vidVaribaleX = solver.AddVariable(null)
let_, vidVaribaleY = solver.AddVariable(null)

//add a row and set it as the goal
let_, vidRow = solver.AddRow(null)
solver.AddGoal(vidRow, 0, true)

letOriginalRosenbrockFunction =
    letf (model:INonlinearModel) (rowVid:int) (values:ValuesByIndex) (newValues:bool) =
        Math.Pow(1. - values.[1], 2.) + 100. * (Math.Pow(values.[2] - (values.[1] * values.[1]), 2.))
    newSystem.Func<INonlinearModel, int, ValuesByIndex, bool, float> (f)

letOriginalRosenbrockGradient =
    letf (model:INonlinearModel) (rowVid:int) (values:ValuesByIndex) (newValues:bool) (gradient:ValuesByIndex) =
        gradient.[1] <- -2. * (1. - values.[1]) - 400. * values.[1] * (values.[2] - (values.[1] * values.[1]))
        gradient.[2] <- 200. * (values.[2] - (values.[1] * values.[1]))
    newSystem.Action<INonlinearModel, int, ValuesByIndex, bool, ValuesByIndex> (f)

solver.FunctionEvaluator <- OriginalRosenbrockFunction
solver.GradientEvaluator <- OriginalRosenbrockGradient
solver.Solve(solverParams);

Console.WriteLine(solver.ToString());

and the output:

> 
Minimize problem
Dimensions = 2
Variable indexes: 
1, 2
Starting point = 0, 0
Solution quality is: LocalOptima
Number of iterations performed: 21
Number of evaluation calls: 27
Finishing point =0.999999997671584, 0.999999995285991
Finishing value =5.74844710320033E-18

val it : unit = ()
>

Implementation for logistic regression

We first get a sample data set online:

module Net = 
    open System.Net

    let fetchUrlSimple (url:string) = 
        let req = WebRequest.Create(url) 
        let response = req.GetResponse()
        use stream = response.GetResponseStream()
        use streamreader = new System.IO.StreamReader(stream)
        streamreader.ReadToEnd()
     
type dataset = 
    { features: float array array; // (instance = float array) array
      mutable labels: int array; // 
    }

    
let iris = 
    let page = Net.fetchUrlSimple @"http://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data"
    let lines = page.Split([|'\n'|], StringSplitOptions.RemoveEmptyEntries)
    let getIrisId (s:string) =     
        if s = "Iris-setosa" then
            2
        elif s = "Iris-versicolor" then
            -1
        else
            1

    let instances = 
        lines
        |> Array.map (fun line -> 
            let s = line.Split ','
            let features = [| "1.0"; s.[0]; s.[1]; s.[2]; s.[3] |] |> Array.map float // add a dumy 1.0 into the fatures 
            let label = getIrisId s.[4]
            features, label
            )
        |> Array.filter (fun (_, label) -> label <> 2)
        
    let F, L = Array.unzip instances
    {
        features = F;
        labels = L;
    }

The logistic regression implementation:

module LogReg = 
    let mutable lambda = 0.1

    let dotProduct (x:float array) (g:ValuesByIndex) =
        let mutable dot = 0.0
        for i=0 to x.Length-1 do
            dot <- dot + x.[i] * g.[i+1]
        dot

    let sigmoid (x:float array) (y:int) (g:ValuesByIndex) = 
        let mutable dot = 0.0
        for i=0 to x.Length-1 do
            dot <- dot + x.[i] * g.[i+1]
        let z = (float y) * dot
        (* code to deal with extream numbers, but with some issues with l(w) calculation
        if z > 30.0 then
            1.0
        elif z < -30.0 then
            0.0
        else
        *)
        1.0 / (1.0 + exp (- z))

    let sigmoid2 (x:float array) (y:int) (g:float array) = 
        let mutable dot = 0.0
        for i=0 to x.Length-1 do
            dot <- dot + x.[i] * g.[i]
        1.0 / (1.0 + exp (- (float y) * dot))

        
    let logregValue(ds:dataset) = 
        let dim = ds.features.[0].Length

        let f (model:INonlinearModel)  (rowVid:int)  (values:ValuesByIndex)  (newValues:bool) = 
            let mutable L = 0.0
            
            for i=1 to dim do
                L <- L + values.[i]*values.[i]
            L <- - (L * lambda / 2.0)
            
            for i=0 to ds.features.Length-1 do
                //L <- L - log (1.0 + exp (- (float ds.labels.[i]) * dotProduct ds.features.[i] values ))
                L <- L + log (sigmoid ds.features.[i] ds.labels.[i] values)

            // printfn "L = %.10f" L
            L
            
        new System.Func<INonlinearModel, int, ValuesByIndex, bool, float> (f)

    let logregGradient(ds:dataset) = 
        // printfn "gradient calculated"
        let dim = ds.features.[0].Length
        let f (model:INonlinearModel) (rowVid:int) (values:ValuesByIndex) (newValues:bool) (gradient:ValuesByIndex) = 
            for j=1 to dim do 
                gradient.[j] <- -lambda * values.[j]
            
            for i=0 to ds.features.Length-1 do
                let coef = (1.0 - sigmoid ds.features.[i] ds.labels.[i] values) * (float ds.labels.[i])
                for j=1 to dim do 
                    gradient.[j] <- gradient.[j] + coef * ds.features.[i].[j-1]

        new System.Action<INonlinearModel, int, ValuesByIndex, bool, ValuesByIndex> (f)



    let makeSolver(ds:dataset) = 

        // set the solver parameters 
        let solverParams = new CompactQuasiNewtonSolverParams()
        //solverParams.IterationLimit <- 11
        //solverParams.Tolerance <- 1e-4
        // the memory parameter m, default 17
        solverParams.IterationsToRemember <- 10

        // solver
        let solver = new CompactQuasiNewtonSolver()

        let dim = ds.features.[0].Length
        // add variables
        for i=1 to dim do 
            solver.AddVariable(null) |> ignore
        
        // add a row and set it as the goal
        let _, vidRow = solver.AddRow(null)
        solver.AddGoal(vidRow, 0, false) |> ignore

        // set funcation evaluator and gradient evaluator
        solver.FunctionEvaluator <- logregValue(ds)
        solver.GradientEvaluator <- logregGradient(ds)
        
        // solve!
        solver.Solve(solverParams) |> ignore
        
        // get the optimized point (w) in the solver
        let w = Array.create dim 0.0
        for i=1 to dim do
            w.[i-1] <- solver.GetValue(i).ToDouble()
        w

Try our logistic regression model on the Iris dataset:

let sol = LogReg.makeSolver(iris)
let nWrongs = ref 0
for i=0 to iris.features.Length-1 do
    let prob = LogReg.sigmoid2 iris.features.[i] iris.labels.[i] sol
    printfn "L = %A sig value = %A label = %A" iris.labels.[i] prob (if prob > 0.5 then "correct" else incr nWrongs; "wrong")

Reference

Tom Mitchell, Generative and Discriminative Classifiers: Naive Bayes and Logistic Regression, new book chapter. This chapter gives a good introduction to logistic regression. For a more theoretic comparison, see: Ng & Jordan: on discriminative vs generative classifiers a comparison of logistic regression and naïve bayes.

Thomas P. Minka, A comparison of numerical optimizers for logistic regression.

It is well written, the first page ready gives a concise introduction to logistic regression. The rest of the paper has a lot of discussion on optimization methods. Tom also has a Matlab package.

Galen Andrew and Jianfeng Gao. Scalable training of L1-regularized log-linear models. In Proceedings of the 24th International Conference on Machine Learning (ICML 2007), pp. 33-40, 2007.

L1-regularized logistic regression with an implementation in C++.

libLBFGS(LBFGS in C): http://www.chokkan.org/software/liblbfgs/

Note

I quite like this post because it combines so many good things together: F#, logistic regression, optimization and playing with sample datasets.

Introducing Sho and using its libraries from F#

2011-02-10T15:38:00.001+08:00

My friend passed me the link for Sho, a data processing and number crunching tool on .Net platform. It is developed in Microsoft Research. I am quite amazed by this tool and like to share a post here.

The frontend of Sho is IronPython. Besides all the standard Python core language and its libraries, Sho also integrates with part of Intel MKL library, which supports linear algebra operations.

There are rumors that there will be a machine learning/data mining library based on Sho in the future.

Key features of Sho

1) IronPython, and its console has some IntelliSense

2) libraries:

* Matrix and linear algebra (Intel MKL©)

both 32bit & 64bit support + multi-threaded

* Statistics (Sho Library)

* Signal processing (Sho Library)

* Optimization (Microsoft Solver Foundation)

3) Interactive visualization

The official site of Sho has some examples and good documents including a reference book.

IronPython?

The language and the open nature make Python a very suitable language for scientific computing. I know that a lot of researchers in data mining area use or occasionally use Python for their research. Three friends of mine are using Python for their research. They are also very good at Matlab. One said to me that he writes his development speed of Python(Scipy/Numpy) and Matlab are about the same. For the running speed, Scipy uses ATLAS + Lapcak, which is about the same as Matlab.

But I think the real power of Python is that your prototype is your deployment. If you ecosystem uses .Net, your number crunching module in IronPython could be ready be used in the system. While for Matlab, deployment is not that easy.

I think IronPython + Sho would be as good as Python+SciPy/Numpy.

I personally don’t quite favor dynamic languages. It is the libraries Sho uses make Sho excellent, F# + Sho libraries rocks too.

Using the Sho API in F#

The following is a sample F# script showing how to use Sho library in F#.

// set the env variable SHODIR to the root of Sho's installation folder
System.Environment.SetEnvironmentVariable("SHODIR", @"C:\Program Files (x86)\Sho 2.0 for .NET 4")


#r @"C:\Program Files (x86)\Sho 2.0 for .NET 4\bin\ShoArray.dll"
#r @"C:\Program Files (x86)\Sho 2.0 for .NET 4\bin\MatrixInterf.dll"

#time

open ShoNS.Array
open System

let rnd = new Random(1)

let vec = new DoubleArray(3)
vec.FillRandom(rnd)

let arr2 = new DoubleArray(3,3)
arr2.FillRandom(rnd)

// big matrix multiplication
let bigArr = new DoubleArray(1200,20000)
bigArr.FillRandom(rnd)
let bigArr2 = bigArr.Transpose()
let bigArrMulti = bigArr * bigArr2
// multi-threaded:
// Real: 00:00:01.745, CPU: 00:00:10.389, GC gen0: 0, gen1: 0, gen2: 0

// svd
let svd = new SVD(arr2)
printfn "U = %A\nD = %A\nV = %A" svd.D svd.U svd.V

// lu
let lu = new LU(arr2)
printfn "L = %A\nU = %A" lu.L lu.U

Notice that Sho linear algebra library is actually licensed from Intel MKL, which is multi-threaded! (Matlab also uses Intel MKL as its basic linear algebra library.)

Basic operations like matrix multiplication is very fast if it uses more cores, as illustrated above.

Note

Fsi.exe runs in 32bit mode. Thus your cannot create large matrices and other objects in the F# interactive. There are workarounds here.

How to quickly learn F#

2011-02-01T11:35:00.001+08:00

I haven’t blogged about F# for sometime now because I am busy with my research work recently. But I nearly use F# everyday to conduct experiments. Also I attracted several of my teammates to get interested in F#, so I write the following guide for my teammates to quickly learn F#.

There are other “learning outline” pages such as the wiki page on Stackoverflow. However the following list provides a path for learning F# and hope to reduce the learning time to be as short as possible.

BTW, it is the Chinese New Year now. Wish everyone a warm Spring Festival!

Step 1

Before reading books/articles, watch the three videos by Luca Bolognese (former F# PM), Luke Hoban (F# PM) and Don Syme (creator of F#).

Step 2

Read the book by Chris Smith (former F# team member, now at Google): Programming F#. After reading the book, quickly scan the following two online books:

· F# Wikibook: http://en.wikibooks.org/wiki/Programming:F_Sharp

· The F# Survival Guide: http://www.ctocorner.com/fsharp/book/default.aspx

They basically cover the same content as Programming F#; reading them gives you a refresh of basic F#.

Step 3

Try to solve the first 50 problems at Project Euler: http://projecteuler.net/ and read the solutions at: http://fsharp-euler.wikispaces.com/ if you have difficulty. You can also search the solutions online; there are many blogs on how to solve Project Euler problems using F#. But before reading the solution, try your best and you can learn F# basics very quickly.

Step 3 alternative

You can skip step 3 (or part of it) if you have experience in building a compiler or interpreter and you are interested in programming languages. The ProgramAsData folder contains the free online book, homeworks (some F# template code is provided) and well-edited slides in the Program as Data course by Prof. Peter Sestoft.

Step 4

Read some chapters of Expert F# 2.0 by Don Syme, and use the book as a reference when you program a bigger/real-world project in F#.

[Optional] Read books like CLR via C# to know more about the .Net platform.

F# Math Provider RE-available

2010-11-22T21:25:00.001+08:00

see:

http://mathprovider.codeplex.com/

WekaSharp: more features

2010-09-19T17:03:00.001+08:00

1. An F# wrapper for Weka. The minimal wrapper in F# for Weka.

2. More features. (this post) This post will contain improvement over the minimal wrapper, e.g. more Dataset processing function, some plot functionality, etc.

3. Tutorial/Usage examples. This post is for end users, who might have no interested in reading the implementation details, but rather knowing how to use this wrapper to perform data mining tasks in .Net.

Here lists a set of features that have been added to WekaSharp.

Parameters for data mining algorithms.

F# features: named parameters and pattern matching.

Data mining algorithm usually have some parameters to set before applying to the data. E.g. the K-means clustering needs the user to provide the number of clusters, the Support Vector Machine classification algorithm need the user to provide the kernel type of the SVM to use, etc.

Also different algorithms have different set of parameters. In Weka, algorithm parameters are encoded as a string (similar to the parameter style of command line tools). However, a first time user would not be able to remember what does “-R 100” mean for a specific algorithm. So we’d like to have a parameter maker for each algorithm. F# fully supports this task with its named parameters.

Let’s see an example. The following code is for making parameters for Logistic Regression:

type LogReg() =
    static member DefaultPara = "-R 1.0E08 -M -1"
    static member MakePara(?maxIter, ?ridge) = 
        let maxIterStr = 
            let m = match maxIter with
                    | Some (v) -> v
                    | None -> -1
            "-M " + m.ToString()
        let ridgeStr = 
            let r = match ridge with
                    | Some(v) -> v
                    | None -> 1.0E08
            "-R " + r.ToString()
        maxIterStr + " " + ridgeStr

and this is an example usage of it:

let logRegPara = LogReg.MakePara (ridge = 1.0)

You can see that besides named parameter functionality, we can also only supply part of the parameters.

Another highlight F# feature in this the Parameter module is pattern matching. Sometimes, the parameters could be quite complicated, e.g. the nearest neighbor search algorithm used in KNN. There are 4 search algorithms, each of them is capable for some of the distance metrics. Discriminate unions and pattern matching make the parameter generating task easy:

    type DistanceFunction = Euclidean | Chebyshev | EditDistance | Manhattan
        with
        static member GetString(distanceFunction) = 
                let d = match distanceFunction with
                        | Some (v) -> v
                        | None -> DistanceFunction.Euclidean
                match d with
                | DistanceFunction.Euclidean -> "weka.core.EuclideanDistance -R first-last"
                | DistanceFunction.Chebyshev -> "weka.core.ChebyshevDistance -R first-last"
                | DistanceFunction.EditDistance -> "weka.core.EditDistance -R first-last"
                | DistanceFunction.Manhattan -> "weka.core.ManhattanDistance -R first-last"

            let distanceFunctionStr = DistanceFunction.GetString distanceFunction 

            let searchAlgorithmStr = 
                let s = match searchAlgorithm with
                        | Some (v) -> v
                        | None -> NNSearchAlgorithm.LinearSearch
                match s with
                | NNSearchAlgorithm.LinearSearch -> 
                    "weka.core.neighboursearch.LinearNNSearch" + " -A " + "\"" + distanceFunctionStr + "\""
                | NNSearchAlgorithm.BallTree ->
                    "weka.core.neighboursearch.BallTree" + " -A " + "\"" + distanceFunctionStr + "\""
                | NNSearchAlgorithm.CoverTree ->
                    "weka.core.neighboursearch.CoverTree" + " -A " + "\"" + distanceFunctionStr + "\""
                | NNSearchAlgorithm.KDTree ->
                    "weka.core.neighboursearch.KDTree" + " -A " + "\"" + distanceFunctionStr + "\""

If you are interested, you can browse the Java code(in Weka’s GUI component) for this task and appreciate how succinct F# gets.

More Dataset manipulating functions

Weka’s dataset class(Instances) is very hard to use. Normally, a dataset is simply a table, with each row representing an instance and each column representing a feature. The dataset class in Weka is too general to be manipulated easily. Weka provides various filters to manipulate datasets. However, these filters are hard to remember their names and applying a filtering needs several lines of code.

I first write a function that transforms a 2 dimensional array into Weka dataset; it can transform the in-memory data more easily into a dataset.

let from2DArray (d:float[,]) (intAsNominal:bool)

The parameter intAsNominal=true means that if a columns are all integers, then treat this column as a nominal column. Otherwise all columns are treated as real ones.

I also wrapper common filters into single functions, that means, applying one filter only uses one line code and they are all put in Dataset module.

The parallel processing

The PSeq module (in F# PowerPack) makes writing computing-bound parallel programs effortless. And in data mining, when it comes to parameter tuning or trying different models on different datasets, it would be fast to run these tasks in parallel on a multi-core machine.

The following code shows two functions performing a bulk of classification tasks:

let evalBulkClassify (tasks:ClaEvalTask seq) = 
    tasks
    |> Seq.map evalClassify
    |> Seq.toList

let evalBulkClassifyParallel (task:ClaEvalTask seq) = 
    task
    |> PSeq.map evalClassify
    |> PSeq.toList

The latter one is the parallel version, which simply change the Seq in the sequential version to PSeq.

Numeric Performance in C, C# and Java

2010-09-05T11:06:00.000+08:00

By occasion, I find the following draft by Prof. Peter Sestoft:

Numeric Performance in C, C# and Java

He implemented several typical numerical programs in C, C# and Java, and tested their running time.

The conclusion in the paper is that managed code is also performant! Read the paper for the details.

Prof. Sestoft is an expert in compiler design and compiler writing. The draft also contains the generated byte code with his annotations. Enjoy!

Just a note:
Although the C code in the paper is native, it does not take the advantage of the modern CPU design. A fully optimized implementation, which is very hard for normal programmers to write, is faster than the native code. This is why Matlab/Numpy's matrix multiplication is faster than your C code. Matlab uses Intel MKL optimized implementation.

WekaSharp: Tutorial for using Weka in F#/.Net

2010-09-03T18:01:00.001+08:00

There are 3 posts in this series:

1. An F# wrapper for Weka. The minimal wrapper in F# for Weka.

2. More features. (not available yet) This post will contain improvement over the minimal wrapper, e.g. more Dataset processing function, some plot functionality, etc.

3. Tutorial/Usage examples. (this post) This post is for end users, who might have no interested in reading the implementation details, but rather knowing how to use this wrapper to perform data mining tasks in .Net.

Installation

The easiest way is to use Visual Studio 2010 with the latest F# PowerPack. Download the release/source code at http://wekasharp.codeplex.com/ and change the dll references to ones at the bin\Release folder.

Visual Studio 2008 should also be able open the project file. Only the parallel functionality cannot be used in 2008 as Parallel Sequence is only added in .Net 4.0.

The support for Mono is not prepared yet. However, it should be easy to figure out.

As all necessary runtimes are also included in the release, so the user does not need to download Weka or IKVM runtimes.

Quick Guide

The user is encourage to look at the script code in Script.fsx and run them to have a general feeling of the wrapper. These examples cover typical usage of WekaSharp.

Reading the other two posts is also very useful for using WekaSharp, and more importantly, for changing the source code of WekaSharp on your own need.

The License

As Weka is GPL, I must make WekaSharp under GPL.

The Efficiency Concern

People may concern how fast is the wrapper compared to the original Weka in Java. I haven’t thoroughly tested this yet. Some casual tests show that they are just the same, at least on Windows platforms. For instance, I run a 3-fold cross validation on the letter dataset using J48 decision trees, both the wrapper and the Weka (run from GUI) use about 20 seconds. It is quite surprising that IKVM’s compiler does so good a job.

As I will show later, the WekaSharp has some parallel constructs, which enable us to utilize multi-core more conveniently.

I’ve provided a module to do profiling. The design is similar to that of Matlab, but better:) You can use Profile.tic() and Profile.toc(“str”) pair to time the execution of F# code. Or You can get multiple timers by using Profile.getWatch().

The Dataset and IO

The Dataset module contains functions to manipulate datasets. Most of the functions are pure, i.e., they don’t change the input dataset and create a new dataset after processing.

Current 4 types of data files are supported: Weka’s ARFF format, LibSvm, Csv and SvmLight. Both loading and saving are supported, i.e. there are 8 functions.

The following code shows how to load an ARFF dataset and set the last column as the label:

let sonar =
   @"D:\temp\datasets-UCI\UCI\sonar.arff"
   |> Dataset.readArff
   |> Dataset.setClassIndexWithLastAttribute

You can also do the two steps in one function:

let iris =
   @"C:\Program Files\Weka-3.6\data\iris.arff"
   |> Dataset.readArffLastAttributeAsLabel

Parameters

Data mining algorithms usually have parameters, sometimes very complicated parameters. F# provides very convenient syntax construct for optional parameters.

E.g. the full set parameters of an SVM has the SVM type, the kernel type, the C value, and the parameters for the kernel, etc. Maybe you just need to set the C value: Parameter.SVM.MakePara(C = c).

Here Parameter is a module for parameters setting. For each data mining algorithm, there are two members:

* DefaultPara. The default parameter string, which is the same as ones used in Weka GUI.

* MakePara. A function to make different parameters.

Here are several more complicated examples of .MakePara method for different algorithms:

Parameter.SVM.MakePara(kernelType = Parameter.SVMKernelType.LinearKernel, C = 10.0)
Parameter.KMeans.MakePara(K=5, maxIter=10, seed=10)

Parameter.KNN.MakePara(distanceFunction = Parameter.Manhattan, K = 3)

and you can also use the IDE support to find which parameters are supported in a data mining algorithm:

Classification and its evaluation

WekaSharp supports most of the common classification algorithms:

type ClassifierType =
   J48 | LogReg | NeuralNet | KNN | NaiveBayes | SVM | LibLinear | LinReg | AdaBoost | Bagging

There are three types of classification tasks:

type ClaEvalTask =
   | CrossValidation of int * Dataset * ClassifierType * parameters
   | RandomSplit of float * Dataset * ClassifierType * parameters
   | TrainTest of Dataset * Dataset * ClassifierType * parameters

To run this task, you need the evalClassify method in Eval module. The following code shows a complete example using J48 as the classifier:

(* playing decision trees on Iris dataset *)
// load the dataset
let iris =
   @"C:\Program Files\Weka-3.6\data\iris.arff"
   |> Dataset.readArff
   |> Dataset.setClassIndexWithLastAttribute

// describe 3 kinds of classification tasks
let j48Tt = TrainTest(iris, iris, ClassifierType.J48, Parameter.J48.DefaultPara)
let j48Cv = CrossValidation(5, iris, ClassifierType.J48, Parameter.J48.DefaultPara)
let j48Rs = RandomSplit(0.7, iris, ClassifierType.J48, Parameter.J48.DefaultPara)

// perform the task and get result
let ttAccuracy = j48Tt |> Eval.evalClassify |> Eval.getAccuracy
let cvAccuracy = j48Cv |> Eval.evalClassify |> Eval.getAccuracy
let rsAccuracy = j48Rs |> Eval.evalClassify |> Eval.getAccuracy

The evalClassify function returns a result object, you can use “.” to that object in the IDE to find out various types of results available. In the above, we use predefined functions to get the accuracy from it.

Clustering and its evaluation

Performing clustering is very similar to classification. You can define two types of clustering:

type CluEvalTask =
   | ClusterWithLabel of Dataset * ClustererType * parameters
   | DefaultCluster of Dataset * ClustererType * parameters

ClusterWithLabel means that you will need to use the label of the dataset to do evaluation. DefaultCluster does not require the dataset has label (actually, it does now allow datasets to have labels either), so the result will only contain the clustering assignments, but not accuracy, etc.

The following code shows a complete clustering example:

let irisLabeled =
   @"C:\Program Files\Weka-3.6\data\iris.arff"
   |> Dataset.readArffLastAttributeAsLabel

let kmeansTask = ClusterWithLabel(irisLabeled, ClustererType.KMeans, Parameter.KMeans.MakePara(K=3))
let emTask = ClusterWithLabel(irisLabeled, ClustererType.EM, Parameter.EM.MakePara(K=3))
let dbscanTask = ClusterWithLabel(irisLabeled, ClustererType.DBScan, Parameter.DBScan.DefaultPara)


let kmeansResult = Eval.evalClustering kmeansTask |> Eval.getClusterSummary
let emResult = Eval.evalClustering emTask |> Eval.getClusterSummary
let dbscanResult = Eval.evalClustering dbscanTask |> Eval.getClusterSummary

Bulk & parallel processing tasks

Sometimes, you need to run multiple tasks. E.g. you need to run the same task multiple times to see the mean and variance of the result, or you need to try different parameters for an algorithm, or you simply have different data mining tasks to run.

The following example shows how to create a bulk of tasks and run them:

// load the data set
let sonar =
   @"D:\temp\datasets-UCI\UCI\sonar.arff"
   |> Dataset.readArff
   |> Dataset.setClassIndexWithLastAttribute

// set different parameters
let Cs = [0.01; 0.1; 1.; 10.; 50.; 100.; 500.; 1000.; 2000.; 5000. ]

// make the tasks with the parameter set
let tasks =
   Cs
   |> List.map (fun c -> Parameter.SVM.MakePara(C = c))
   |> List.map (fun p -> CrossValidation(3, sonar, ClassifierType.SVM, p))
  
Profile.tic()
// the accuracy result
let results =
   tasks
   |> Eval.evalBulkClassify
   |> List.map Eval.getAccuracy
Profile.toc("sequential time: ")

Here I created different SVM tasks for different C values, run them and get the accuracy as a list.

F# provides very easy syntax to perform multiple tasks at the same time. Thus I provide evalBulkClassifyParallel method:

Profile.tic()
let resultsParallel =
   tasks
   |> Eval.evalBulkClassifyParallel
   |> List.map Eval.getAccuracy
Profile.toc("parallel (PSeq) time: ")

// sequential time: : 9767.804800 ms
// parallel (PSeq) time: : 6154.715500 ms

As the profiling shows, on a two-core machine, parallel executing the tasks does boost the speed.

Plotting

is not finished, but it is still quite usable. To plot the different accuracies for different Cs in the SVM, we can use:

// do the plot
lc.column(y = results, xname = "differnet C", yname = "Accuracy", title = "SVM on iris",
   isValueShownAsLabel = true ) |> display

Making Datasets from Memory

All the above examples use data from data files. In practice, we might want to convert the data in memory, e.g. in an array into a Weka dataset. The following examples shows how to create dataset from F# arrays:

(* create dataset from F# arrays *)

// make the data array
let data = [| 0.; 0.;
             1.; 1.;
             0.; 1.;
             1.; 0.; |]
let xorArray = Array2D.init 4 2 (fun i j -> data.[i*2 + j])

// make weka dataset from array
let xor0 = Dataset.from2DArray xorArray false

// add labels
let xor = xor0 |> Dataset.addClassLabels ["T"; "T"; "F"; "F"]

// make svm tasks

let rbfTask = TrainTest(xor, xor, ClassifierType.SVM, Parameter.SVM.DefaultPara)
let linearTask = TrainTest(xor, xor, ClassifierType.SVM, Parameter.SVM.MakePara(kernelType = Parameter.SVMKernelType.LinearKernel, C = 10.0) )

// rbf svm gets 100% accuracy
let rbfAccuracy = rbfTask |> Eval.evalClassify |> Eval.getAccuracy
// linear svm does not work on XOR data set
let linearAccuracy = linearTask |> Eval.evalClassify |> Eval.getAccuracy

Conclusion

As a user of WekaSharp, I already benefit from it as I can process the data and run data mining algorithms over it both in F#. Originally, I needed to write my data into ARFF file, call Weka, parse the output to get the result. And the F# solution is far more declarative.

Another benefit is that we can write extensions to Weka in F#!

Enjoy!

WekaSharp: An F# wrapper for Weka

2010-08-23T20:00:00.001+08:00

There are 3 posts in this series:

1. An F# wrapper for Weka. (this post) The minimal wrapper in F# for Weka.

2. More features. This post will contain improvement over the minimal wrapper, e.g. more Dataset processing function, some plot functionality, etc.

3. Tutorial/Usage examples. This post is for end users, who might have no interested in reading the implementation details, but rather knowing how to use this wrapper to perform data mining tasks in .Net.

WekaSharp is available at CodePlex:http://wekasharp.codeplex.com/.

Weka is one of the most widely used data mining software package. It contains over 50 data mining algorithms, a good GUI support and well written documents. Because of its good features, a lot data mining courses use it as an illustrative software.

Its GPL license also allows it to be used freely in academic.

However, Weka GUI enables us to perform data loading, filtering and classifier building in an interactive way, sometimes programming is still needed. For instance, GUI has poor support for parameter selection and processing 100 data sets. In these tasks, scripting is a better way. However, Java language is not declarative and it does not support interactive shell either.

F# well supports these two features, also data processing in F# in general is far better than Java. So F# and Weka are good combination.

Compile Weka library to a .Net dll

As everybody knows, Weka is written in Java. To use it within .Net framework seamlessly, we need to compile it into a .Net assembly. Fortunately, IKVM project is made for such purpose.

There is an online tutorial for compiling weka.jar to weka.dll. However, this tutorial does not mention how to use external jar libraries in Weka, e.g. libsvm.jar. Because Java’ class loading mechanism is different from .Net’s. We cannot simply compile different jars into separate dlls. The simplest way is to compile weka.jar and libsvm.jar to a single dll:

>ikvmc -sharedclassloader -target:library weka.jar libsvm.jar

Notice that there is a WekaDotNet project at sourceforge. The project is compiled using Visual J#, which is already deprecated in Visual Studio family. Using IKVM to compile jar into .Net is an easier and more stable way.

The overall structure of the wrapper

Based on the Weka class hierarchy, I’ve defined several separate modules for different kinds of functionality:

* Common – contains all common types over the whole wrapper. E.g. the union type for various classifiers.

* Dataset – contains IO routines to read/save various datasets and functions (e.g. split) to manipulate/preprocess datasets.

* Parameter – contains functions to make parameters for the data mining algorithms, all the default parameters are also provided.

* Classify – contains standard classifiers.

* Cluster – contains standard clustering algorithms.

* Eval – contains functions to perform different kinds of classification/clustering tasks and output their result (e.g. accuracy/F1 measure) out.

In the following sections, I will present each module in details.

Common module

In this module, I try to write all types, including classifier types, clustering algorithm types and some evaluation-task types. This module is marked as AutoOpen, i.e., the types in the modules are easily accessible from other modules.

[<AutoOpen>]
module Common =
let NYI() = failwith "Not yet implemented!"

type Dataset = core.Instances

type parameters = string

type DatafileType =
   Arff | LibSVM | Csv | Svmlight

type ClassifierType =
   J48 | LogReg | NeuralNet | KNN | NaiveBayes | SVM | LibLinear

type ClaEvalTask =
   | CrossValidation of int * Dataset * ClassifierType * parameters
   | RandomSplit of float * Dataset * ClassifierType * parameters
   | TrainTest of Dataset * Dataset * ClassifierType * parameters

type ClustererType =
   KMeans | EM | DBScan

I’ve redefined Weka Instances class as Dataset, and use a string for parameters.

Dataset module

This module includes functions to read/save 4 kinds of data files: Weka’s ARFF, libsvm, Csv and SvmLight. It also has a function randomSplit to randomly split a dataset with a given ratio.

This module currently is far from complete. For instance, we can only load datasets from 4 kinds of disk files, and no way to build a dataset using in-memory data, e.g. an array of vectors. We also only provide one data preprocessing step – random split, many other common preprocessing steps are needed.

I’d like to use another post to finish these. In this post, delivering the whole working wrapper (although incomplete) is more important.

Here is some of the implementation:

    let readDataFile (ftype: DatafileType) (fname:string) =
   let loader =
       match ftype with
           | Arff -> new converters.ArffLoader() :> AbstractFileLoader
           | LibSVM -> new converters.LibSVMLoader() :> AbstractFileLoader
           | Csv -> new converters.CSVLoader() :> AbstractFileLoader
           | Svmlight -> new converters.SVMLightLoader() :> AbstractFileLoader
   try
       loader.setFile(new java.io.File(fname))
       loader.getDataSet()
   with
       | _ -> failwith "dataset loading error"

its input is a datafiletype and the file name. Based on this function, 4 concrete functions are defined:

let readArff = readDataFile DatafileType.Arff

    let readLibsvm = readDataFile DatafileType.LibSVM
    let readCsv = readDataFile DatafileType.Csv
    let readSvmlight = readDataFile DatafileType.Svmlight

The data saving functions have similar implementation.

Parameter module

This module contains functions that create parameters for different data mining algorithms. As Weka uses a space-separated string as parameters, it would be unclear to the first-time user that what does “-R 1.0 –M -1” mean for a logistic regression.

I have provided a default parameter string for each algorithm.

Classify module

In this module, I have wrapped most commonly used classification algorithms: C4.5 decision tree (named J48 in Weka), NaiveBayes, Logistic Regression, Multilayer perception neural net, Support Vector Machines (SVM).

The getClassifier function accepts a classifier type and the an option string, and returns a classifier instance:

let getClassifier (ctype:ClassifierType) (option:string)  =
let classifier =
   match ctype with
       | J48 -> new classifiers.trees.J48() :> classifiers.Classifier
       | NaiveBayes -> new classifiers.bayes.NaiveBayes() :> classifiers.Classifier
       //| KNN ->new classifiers.lazy.IBk() :> weka.classifiers.Classifier
       | LogReg -> new classifiers.functions.Logistic() :> classifiers.Classifier
       | NeuralNet -> new classifiers.functions.MultilayerPerceptron() :> classifiers.Classifier
       | SVM -> new classifiers.functions.LibSVM() :> classifiers.Classifier
       | LibLinear -> new classifiers.functions.LibLINEAR() :> classifiers.Classifier
       | _ -> failwith "not supported"

classifier.setOptions(core.Utils.splitOptions option)
classifier

and buildClassifier takes one more parameter – the dataset to build a trained classifier:

let buildClassifier (ctype:ClassifierType) (option:string) (ds:Dataset) =
checkDataset ds
let classifier = getClassifier ctype option
classifier.buildClassifier(ds)
classifier

Handy functions like getJ48, getSVM are also defined:

let getJ48 (option:string) = getClassifier ClassifierType.J48 option
let getSVM (option:string) = getClassifier ClassifierType.SVM option

There are several issues:

1. KNN (Weka class: IBk) is not supported as its namespace in Weka is core.classifiers.lazy. But lazy is a keyword in F#, thus cannot appear in a namespace.

2. The IO wrapper for classifiers are not provided yet.

Cluster module

Similar to Classify module, there are getClusterer and buildClusterer in Cluster module. There are also shortcuts to concrete clustering algorithms such as getKmeans and buildKmeans.

Eval module

This module contains evaluation methods for classification and clustering. For classification, I’ve defined a task discrete union:
type ClaEvalTask =

    | CrossValidation of int * Dataset * ClassifierType * parameters
| RandomSplit of float * Dataset * ClassifierType * parameters
| TrainTest of Dataset * Dataset * ClassifierType * parameters

classifyEval function is to do such a task:

let rec classifyEval (task:ClaEvalTask) =
let rnd = new java.util.Random(System.DateTime.Now.Ticks)
match task with
| TrainTest (dsTrain, dsTest, ctype, para) ->
   Classify.checkDataset dsTrain
   Classify.checkDataset dsTest
   let eval = new classifiers.Evaluation(dsTrain)
   let cl = Classify.buildClassifier ctype para dsTrain
   eval.evaluateModel(cl, dsTest) |> ignore
   eval
| CrossValidation(cv, ds, ctype, para) ->
   Classify.checkDataset ds
   let eval = new classifiers.Evaluation(ds)
   let cl = Classify.getClassifier ctype para
   eval.crossValidateModel(cl, ds, cv, rnd, Array.empty)
   eval
| RandomSplit(ratio, ds, ctype, para) ->
   Classify.checkDataset ds
   let train, test = Dataset.randomSplit ratio ds
   classifyEval (TrainTest(train, test, ctype, para))

the result of this function is an evaluation result object. We can use its properties to get various evaluation metrics, or use the shortcut functions in the Eval module:

let getAccuracy (eval:classifiers.Evaluation) = eval.pctCorrect()

let getPrecison (eval:classifiers.Evaluation) (classIdx:int) = eval.precision(classIdx)
let getRecall (eval:classifiers.Evaluation) (classIdx:int) = eval.recall(classIdx)
let getF1 (eval:classifiers.Evaluation) (classIdx:int) = eval.fMeasure(classIdx)
let getAUC (eval:classifiers.Evaluation) (classIdx:int) = eval.areaUnderROC(classIdx)
let getClassifySummary (eval:classifiers.Evaluation) = eval.toSummaryString()

Conclusion

Two features of the wrapper are

1) a declarative wrapper to Weka. See the following code to see how declarative the wrapper is:

// load the dataset

let dataset =
@"C:\Program Files\Weka-3.6\data\iris.arff"
|> Dataset.readArff
|> Dataset.setClassIndexWithLastAttribute

// describe the cross validation task using SVM as the classifier
let classifyTask = CrossValidation(5, dataset, ClassifierType.SVM, Parameter.SVM.defPara)

// perform the cross validation task and then get the average accuracy
let cvAccuracy =
classifyTask
|> Eval.classifyEval
|> Eval.getAccuracy

The code is understandable even to those who do not know F#.

2) interactively perform data mining algorithms.

F# interactive is well integrated in Visual Studio. It is very stable and convenient to use. The powerful VS IDE provides good Intellisense for F#.

Reading F# Projects, Part III: The F# math providers.

2010-08-14T21:25:00.001+08:00

We’ve covered the usage of Math providers in the matrix series. In this post, the whole source code is decomposed.

As I find it time consuming to covert One Note notes into a blog format. Please see the attached pdf file.

PDF.

Reading F# Projects, Part II: F# Set

2010-08-14T20:49:00.001+08:00

The design patterns in a data structure project like F# Set is simple. However, understanding it really needs time as data structures involve heavy generic programming and the type system in F# (or any other static FP language) is very complex.

In this post, we will explore two aspects of F# Set:

1. The general implementation structure for a generic data structure. After knowing the pattern, we could write other data structures in F# and these data structures are generic enough.

2. Know something about generics in F# and type inference. We will come to this part in later readings.

In F# source code, the immutable set implementation (set.fs) is a good one for us to start for data structures. Why not more fundamental data structures like F# sequences? Because those ones are more involved in .Net platform and F# Set is nearly pure F# code. Also because:

I love trees!

F# Set and Map are implemented using balanced binary search tress, to be exact, AVL tress. Balanced BSTs support insert, delete and find operations all in O(log n). So it is widely used for implementing dynamic sets.

In general, trees are the most beautiful data structures and have so many variants. They are also efficient, e.g. R tree and friends are the fundamental data structures used in every database. Quadtrees are used for computational geometry in theory and spatial database in practice.

The F# Set is implemented by an immutable AVL tree. Here immutable means that every node, once created, its values could be be changed. In CS literature, researchers use “persistent” to name immutable data structures. If you do 100 insertions to an empty tree, you can have 101 version of trees, where common nodes are shared to reduce the space cost. It is proved that in such a tree, inserting operation still costs O(log n) time, and O(log n) space for duplicating some nodes in the tree (usually along the insertion path).

As this post is not intended to be a tutorial on data structures, I won’t go into details for the tree implementation. We mainly study how an industry-strength data structure is implemented in F#.

The main structure

set.fs contains about 840 lines of code. The code skeleton is shown below:

namespace Microsoft.FSharp.Collections
   type SetTree<'T> when 'T : comparison =
       | SetEmpty                                          // height = 0  
       | SetNode of 'T * SetTree<'T> *  SetTree<'T> * int    // height = int
       | SetOne  of 'T                                     // height = 1  
   module internal SetTree =
       // functions:
       // indert, delete, find, etc.
       // split, rebalance, etc.
      
       type SetIterator<'T> when 'T : comparison  =
           { mutable stack: SetTree<'T> list;  // invariant: always collapseLHS result
             mutable started : bool           // true when MoveNext has been called  
           }
       type Set<'T when 'T : comparison >(comparer:IComparer<'T>, tree: SetTree<'T>) =
           //member functions:
           //  Add, Remove, Contains, etc.
           //operators:
           //  - + Intersection, etc.
           //overrides:
           //  GetHashCode, Equals, ToString
           //implement interfaces:
           //  IComparable
           //  ICollection<'T>
           //  IEnumerable<'T>
           //  IEnumerable
   module Set =
       // set module functions:
       // add, isEmpty, etc.

Under Microsoft.FSharp.Collections namespace, there are three things: SetTree<’T> type, Set module and SetTree module.

SetTree is a type constructor, the parameter is a comparable type ‘T.

Most of the code is in SetTree internal module:

the first part contains the code to manipulate the tree structure, e.g. inserting a node, deleting a node, etc.

Type SetIterator<’T> and related functions are for implementing IEnumerable interface.
Set<’T> (comparer, tree) is the implementation for Set.

Set module contains the common function to manipulate a set.

In a summary, A container, like Set, usually follows this pattern:

1. An implementation of the underlying data structure, here it is AVL tree.

2. A module for manipulating the data structures; members functions also can do this. Usually we keep functions in both places.

3. Implement the IEnumerable and IEnumerable<’T> interface, so that the container/data structure could be treated as a sequence.

These are kind of rules for implementing new containers in F#, e.g. F# Map.

Detail 1. The generic

It is meaningless if the Set implementation only supports integers, or floats. We want it to be generic, i.e., support any type. Since the underlying implementation is based on trees, so the base type should also support comparison.

The other part of generic, which is more subtle, is the comparer. There are situations, we have different string sets. However, the rules to order them are not the same, i.e., strings maybe compared using different rules. So we should add another parameter into consideration. (Refer to Functor in OCaml if you want to see a more elegant more to deal with this situation.)

Although F# Set currently does not support different comparers. Its internal implementation keeps comparer in mind. (Because F# supports OOP, the different comparers problem could be partially solved by inherit the original element type. Anyway, this is non-elegant.)

We know that once a Set has some concrete values in it, the type is clear. Here are two ways to initialize a Set:

static let empty : Set<'T> =

    let comparer = LanguagePrimitives.FastGenericComparer<'T>
   new Set<'T>(comparer, SetEmpty)

new (elements : seq<'T>) =
   let comparer = LanguagePrimitives.FastGenericComparer<'T>
   new Set<_>(comparer,SetTree.ofSeq comparer elements)

We can see that they get the comparer for type ‘T and use that comparer to pass into the constructor of type Set<’T>.

Detail 2. The comparer

The above code get the comparer for type ‘T using FastGenericComparer, which is in a big file prim-types.fs. The primitive type module should be discussed separately and thoroughly. There are quite a few ad hoc tricks in the implementation.

Here we just need to know that it looks up a table to find a comparer and if the type is in the primitive types, then there are optimizations for them.

Detail 3. From the empty to concrete

F# allows you to write [], Set.Empty. These empties seem to be dynamic. Actually no. You cannot write sth like:

let empties = Array.create 100 Set.empty

the compiler gives the following error:

error FS0030: Value restriction. The value 'empties' has been inferred to have generic type

    val empties : Set<'_a> [] when '_a : comparison

Because the compiler cannot generalize these empties into a single type. It is safe to write

Set.empty + Set([1;2;3])

Because the types of the two sets could be generalized.

Same situation applies to the zeros in different number types.

Detail 4. The Iterator

interface IEnumerable<'T> with
   member s.GetEnumerator() = SetTree.mkIEnumerator s.Tree

interface IEnumerable with
   override s.GetEnumerator() = (SetTree.mkIEnumerator s.Tree :> IEnumerator)

Both of the IEnumerable interfaces are implemented. Notice that difference between the two: there is a upcast in the type. Let’s see the actual implementation:

let mkIterator s = { stack = collapseLHS [s]; started = false }

let mkIEnumerator s =

    let i = ref (mkIterator s)
   { new IEnumerator<_> with
         member x.Current = current !i
     interface IEnumerator with
         member x.Current = box (current !i)
         member x.MoveNext() = moveNext !i
         member x.Reset() = i :=  mkIterator s
     interface System.IDisposable with
         member x.Dispose() = () }

Exercise

An exercise is to read map.fs in F# core. The implementation of F# Map is very similar to F# Set.

Reading F# Projects, Part I: The Common Knowledge

2010-08-14T20:48:00.001+08:00

During the programming of my data mining library, I occasionally refer to the design of other libraries. E.g. data mining toolkits WEKA and MALLET are well designed. I also refer to the design of F# libraries. The examples on F# books mainly tell us how to program on small level, i.e. how to manipulate list or sequence and some language features. However, designing in F# is often covered in a general way.

Thus, it is beneficial for new F# users to read well designed F# libraries. In this series, the design of the following source code will be presented: (the list is tentative to change.)

* some parts of F# Core and PowerPack.

* F# math providers -- service design pattern and PInvoke in F#.

* TrueSkill – how a specific data mining model is implemented.

* FParsec -- the best example of computation expression/monad.

Before going to any concrete projects, we review the common issues/tips/concepts that may be involved in the process of the project design in F#.

I will add more topics here as we go through different F# programs.

The interface/signature file .fsi

In Java or C#, the implementation and interface are put in one. In C/C++, the interface (.h) and the implementation (.c/c++) are separated. Either design is good and has advantages.

For F#, I think the later is better. In .fsi, all the public functions/types are commented carefully. These commented can be used to generated documents for the program, or used by the IDE Intellisense.

The actual implementation is in the .fs files. In this file, comments are in places where actually necessary for the program who implements it, e.g. a formula, or a refer to a page number in a book. Because the program who writes the function definitely know what every parameter mean. If this kind of comments are put into the implementation, then the programs look less succinct.

So usually, the programmer doing the implementation write the implementation in .fs. When the implementation is stable or ready to ship out, the corresponding .fsi interface file is generated and well commented manually.

Here is the MSDN page on signature file.

Namespaces and Modules

The latest F# version requires that an F# source file (.fs) starts with a namespace, that means all the following code is under this namespace(An exception: there could be multiple namespaces in a source file.)

Modules are similar to namespaces. You can view that a group of functions are put into a namespace.

The usage difference between namespaces and modules, in my option, is that namespace is broader concept. E.g. Microsoft.FSharp.Collections names contains a lot of standard data structures, where the manipulating functions for each data structure are put into a module.

You can also view modules as a class containing only static members.

Classes and Interfaces

F# interfaces are just like interfaces in C#/Java. Interface is a quite universal concept occurring different programming paradigms. Even in the pure FP Haskell, type classes are similar to interfaces.

F# classes are different from classes in OO languages. First, F# encourages an immutable programming fashion for classes: you have only one main constructor, once the object is constructed the values remain unchanged. You can also have mutable member fields in a class. The constraints on F# classes are for safer programs, although sometimes causes some inconvenience.

Besides classes, F# also has Records, Enums and Discrete Unions.

Extension to existing classes and modules

In F#, you can extend an existing .net class by using “type .. with ..” construct:

type System.Net.WebRequest with
   member this.AsyncGetResponse() : Async<System.Net.WebResponse>=
       async {… the implementation }

(from F# source code: control.fs)

In this piece of code, System.Net.WebRequest is an existing class in .Net. We add a member function AsyncGetResponse into this class. This is like lightweight inheritance by saving from creating a new class.

The attributes

In .Net, attributes associate declarative information with the code. Usually they occur in production/formal code. Here is a tutorial.

Here I list some commonly used attributes in F#:

[<CompiledName("xxx")>]
[<AutoOpen>]
[<RequireQualifiedAccess>]

To see a full list, and the exact meaning. No resource is better than the F# Specification, Chapter 16. Special Attributes and Types.

F# Async workflow application: a Flickr crawler

2010-07-10T23:58:00.001+08:00

In data mining, much time is spent on collecting raw data and preprocessing data. Quite a few data mining research tasks require to download data from Internet, e.g. Wikipedia articles from Wikipedia, photos from Flickr, Google search results, reviews from Amazon, etc.

Usually general purpose crawlers, such as wget, are not sufficiently powerful and specialized in downloading data from Internet. Writing a crawler on our own is often required.

Recently I am doing an image related research project, in which I need to download a lot of tagged images from Flickr. I am aware that there is a Flickr downloadr, which uses Flickr API to download images. However 1) it only downloads licensed photos and 2) it cannot download the tags of a photo. Thus I decided to write one myself.

The input of the is a tag query, e.g. “dog”, # of photos to download and the disk folder to store the downloaded images.

Because the number of photos is quite big, so downloading them in parallel is critical. In .Net, there are several ways to do parallel computing. For IO intensive tasks, Async workflow is the best.

Tutorials for Async workflow

As Async workflow is one of the key features of F#, there are quite a few tutorials online for F# Async programming. Providing one more in this blog would be repetious.

Luke Hoban’s PDC 2009 talk, F# for Parallel and Asynchronous Programming, is very good for beginners.

Don Syme wrote 3 articles in a series. Highly recommended for experienced F# users!

The Flickr crawler

To write a crawler for a web site like Flickr, we need to 1) design the downloading strategy and 2) analyze the structures of Flickr.

My strategy is to use the search query to search images with some specific tags and from the result page(as shown below), the url of each image is extracted, from which the image and its tags are then crawled.

Figure 1. Flickr search result page.

Figure 2. Flickr image page with an image and its tags.

So first we need a function to download a web page asynchronously :

 let fetchUrl (url:string) = 
    async {
        try 
            let req = WebRequest.Create(url) :?> HttpWebRequest
            req.UserAgent <- "Mozilla/4.0 (compatible; MSIE 6.0; Windows NT 5.1; SV1; .NET CLR 1.1.4322)";
            req.Method <- "GET";
            req.AllowAutoRedirect <- true;
            req.MaximumAutomaticRedirections <- 4;
            let! response1 = req.AsyncGetResponse()
            let response = response1 :?> HttpWebResponse
            use stream = response.GetResponseStream()
            use streamreader = new System.IO.StreamReader(stream)
            return! streamreader.AsyncReadToEnd() // .ReadToEnd()
        with
            _ -> return "" // if there's any exception, just return an empty string
    }
 

fetchUrl pretends to be a Mozilla browser and can do some redirections if the url is slightly invalid. The current exception handling is very easy – just return the empty string for the web page. Notice that the return type of the function is Async<string>, thus it cannot be used to download images as images are of binary format, not text.

So the next task is to write a function to download images:

 let getImage (imageUrl:string) =
    async {
        try 
            let req = WebRequest.Create(imageUrl) :?> HttpWebRequest
            req.UserAgent <- "Mozilla/4.0 (compatible; MSIE 6.0; Windows NT 5.1; SV1; .NET CLR 1.1.4322)";
            req.Method <- "GET";
            req.AllowAutoRedirect <- true;
            req.MaximumAutomaticRedirections <- 4;
            let! response1 = req.AsyncGetResponse()
            let response = response1 :?> HttpWebResponse
            use stream = response.GetResponseStream()
            let ms = new MemoryStream()
            let bytesRead = ref 1
            let buffer = Array.create 0x1000 0uy
            while !bytesRead > 0 do
                bytesRead := stream.Read(buffer, 0, buffer.Length)
                ms.Write(buffer, 0, !bytesRead)
            return ms.ToArray();
 
        with
            _ -> return Array.create 0 0uy // if there's any exception, just return an empty image
    }
 

Next we write the code to get the url of an image and its tags from an image page(see Figure 2):

 let getBetween (page:string) (head:string) = 
    let len = head.Length
    let idx = page.IndexOf(head)
    let idx2 = page.IndexOf('"', idx+len)
    let between = page.Substring(idx+len, idx2 - idx - len)
    between
 
 
let getImageUrlAndTags (page:string) = 
    let header = "class=\"photoImgDiv\">"
    let idx = page.IndexOf(header)
    let url = getBetween (page.Substring(idx)) "<img src=\""
 
    let header2 = "<meta name=\"keywords\" content=\""
    let tagStr = getBetween page header2
 
    let s = tagStr.Split([|','|], System.StringSplitOptions.RemoveEmptyEntries)
    let tags = s |> Array.map (fun t -> t.Trim())
    url, tags
 

Finally, write a function to work through every search result page, parse the result page and download the images in that result page:

 let getImagesWithTag (tag:string) (pages:int) = 
    let rooturl = @"http://www.flickr.com/search/?q="+tag+"&m=tags&s=int"
    seq {
        for i=1 to pages do 
            let url = rooturl + "&page=" + i.ToString() 
            printfn "url = %s" url
            let page = fetchUrl url |> Async.RunSynchronously
            let imageUrls = getImageUrls page
            let getName (iurl:string) = 
                let s = iurl.Split '/'
                s.[s.Length-1]
            
            (* images in every search page *) 
            let images = 
                imageUrls 
                |> Seq.map (fun url -> fetchUrl url) 
                |> Async.Parallel
                |> Async.RunSynchronously
                |> Seq.map (fun page -> 
                    async {
                        let iurl, tags = getImageUrlAndTags page
                        let icontent = getImage iurl |> Async.RunSynchronously
                        let iname = getName iurl
                        return iname, icontent, tags
                    })
                |> Async.Parallel
                |> Async.RunSynchronously
            yield! images
        }
 

with a driver function to write all the images into hard disk:

 let downloadImagesWithTag (tag:string) (pages:int) (folder:string) = 
    let images = getImagesWithTag tag pages
    images
    |> Seq.iter (fun (name, content, tags) -> 
        let fname = folder + name
        File.WriteAllBytes(fname, content)
        File.WriteAllLines(fname + ".tag", tags)
        )
 

We’ve done! A Flickr image crawler in only about 120 lines of code. Let’s download some images!

downloadImagesWithTag "sheep" 5 @"D:\WORK\ImageData\flickr\sheep\"

It only costs less than 5 minutes to download 300 sheep pictures from Flickr.

Discussions

1. One of the strengths of F# Async programming is its ease for exception handling. In this example, the exception is handled immediately, we can also propagate the exception into an upper level function and handle it there. However, that would require more thinking….

2. I only parallelly download images in one search page. The program could be modified to parallelly process multiple search result pages, which is done sequentially now. If done in this way, we can see that, we can build a hierarchical parallel program: 1) at the first level, multiple search result pages are processed parallelly and 2) at the second level, images in a research result pages are downloaded parallelly.

A note on efficiency in F#, part I: small things first

2010-06-04T22:42:00.001+08:00

In this blog, I often talk about efficiency aside data mining as efficiency really matters in this field. Previously, we’ve talked about how to use PInvoke and C++/CLI to use C++ code to boost the performance. In this post, I’d like to start a series on F# itself. As a functional programming language, F# programs are kind of high-level, thus harder to reason its performance than imperative programs. Not the big-O thing, it is the constant. This series contains little tips about performance.

Note that these tips are just facts to let you know. There is no performance guidance that one should follow. There is always a treadoff, high level functional constructs usually are slower than their pure imperative equivalents. However, they are safer and have better abstraction. It usually depends when you really need to do some optimization. E.g. if you could guarantee that your program runs less than 0.5 second, then optimizing it would usually become less meaningful, unless you have another program calling this program a million times.

Also remember that

"We should forget about small efficiencies, say about 97% of the time: premature optimization is the root of all evil" (by Donald Knuth)

1. Stream IO and Seq module

Think about that we need to process a large text file, say, 1TB. Then it is infeasible to load the file into memory. The approach is to treat the text file as a stream and scan this stream to get useful information.

 open System.IO
let readLinesSeq = File.ReadLines
 

.Net 4.0 provides ReadLines, the type of its return value is seq<string>.

In lower versions of .Net, one needs to implement this function:

 let readLines filePath = seq { 
    use sr = new StreamReader (filePath) 
    while not sr.EndOfStream do 
        yield sr.ReadLine () 
}
 

Once we have this function, we can process it as as a sequence.

While prototyping a model, we may make mistakes, and we also know that scanning the whole file would cost a lot of time. If we find an error after scanning the whole scanning, that would waste a lot of time. So we would like to take the first few lines of the file as a small test set for the debugging purpose. By treating a file as a sequence, we could simply use

 let lines100 = "filename" |> readLines |> Seq.take 100
 

to get a small test set.

Here is an quite old post in 2006. At that time, there seemed to be no Seq module. The idea is the same.

2. Loop, Seq module and the tail-recursive version of “while”

One of my first F# programs is a prime number test function. I first used an ugly for loop with mutable variables:

 let isPrime1 n = 
    let m = int(sqrt (float n))
    let mutable flag = true
    for i=2 to m do 
        if n % i = 0 then 
            flag <- false
    flag
 

You can see that I clearly didn’t know how to early stop/break the loop when flag is assigned to false. Because I could not find break, so I crazily use exception to break out the loop:-)

 exception Break 
let isPrime1' n = 
    let m = int(sqrt (float n))
    let mutable flag = true
    try
        for i=2 to m do 
            if n % i = 0 then 
                flag <- false
                raise Break
    with
        _ -> ()
    flag
 

And based on my FP experience on Scheme, I immediately wrote a tail-recursive function, which behaves like a “while” loop:

 let isPrime2 n = 
    let m = int(sqrt (float n))
    let rec loop div = 
        if div > m then true
        elif n % div = 0 then false
        else loop (div+1)
    loop 2
 

I can also explicitly use a while loop with a mutable flag:

 let isPrime3 n = 
    let m = int(sqrt (float n))
    let mutable flag = true
    let mutable div = 2
    while flag && div <= m do 
        if n % div = 0 then
            flag <- false
        div <- div + 1
    flag
 

And I finally find the elegant way using Seq:

 let isPrime n = {2..int(sqrt (float(n)))} |> Seq.forall (fun i->n%i<>0)
 

However, the last version would be slower than others except for the first one. Because a sequence is actually an object implements IEnumerable interface. If the object currently points to 10, to get 11, it needs to call a Next function to get it. There is definitely some overhead here.

So pay attention when you use Seq to do the job of a for loop or while loop when the loop is executed for a huge number of times.

3. List and Array

F# lists are single linked lists. Thus, Lists use more memory as it needs to keep a next pointer in every node. Lists are also slow. Let’s do a test:

 let n = 10000000
let a = List.init n (fun i-> i)
let b = Array.init n (fun i-> i)
 

creating list a uses about 2.2 seconds and about 200MB memory, while creating array b only uses 0.12 second and about 40MB memory, which equals our estimate n * 4 bytes/int ~= 40MB.

4. Access lists, arrays and sequences

We know that the three modules List, Array and Seq have similar functions. E.g. all of them have the function map, which applies a function to the existing list, array or sequence to generate a new one. However, the three maps are implemented differently, with array most efficiently, list second, the most general seq the last:

 // form F# source code: local.fs 
// note: the actuall implementation for list is in local.fs, not in list.fs
let map f x = 
    match x with
    | [] -> []
    | [h] -> [f h]
    | (h::t) -> 
        let cons = freshConsNoTail (f h)
        mapToFreshConsTail cons f t
        cons
// form F# source code: array.fs 
let map (f: 'T -> 'U) (array : 'T[]) : 'U[]=
    checkNonNull "array" array
    let inputLength = array.Length
    let result = Microsoft.FSharp.Primitives.Basics.Array.zeroCreateUnchecked inputLength
    Parallel.For(0, inputLength, fun i ->
        result.[i] <- f array.[i]) |> ignore
    result
// form F# source code: seq.fs 
let map f (e : IEnumerator<_>) : IEnumerator<_>=
    upcast 
        { new MapEnumerator<_>() with
            member this.DoMoveNext (curr : byref<_>) =
                if e.MoveNext() then
                    curr <- (f e.Current)
                    true
                else
                    false
            member this.Dispose() = e.Dispose()
        }
 

We can see that the Array.map is actually implemented using ParallelFor. The most general implementation is Seq.map based one the IEnumerable interface in .Net, which is relatively inefficient compared to a for loop.

However, by using Array.map or List.map, we create another array or list. While Seq.map only creates an auxiliary sequence, which costs nearly no memory!

So usually, I use Seq.map to map arrays or lists to keep memory usage small. When performance really matters, I use array and for loop, i.e. I write imperative code. In my experience, writing imperative code in F# sometimes still introduces some runtime errors, while this is rare in pure functional code.

When I only need to iterate an array, I will directly use Array.iter since it is efficient, and use no extra memory.

All in all: performance costs. Do some quick estimation before you write imperative code. Also know the low-level implementation of your daily-use libraries is critical for you to estimate the performance.

5. Access F# matrices and vectors

The Vector<’T> type is implemented using 1-D array and the Matrix<’T> type is implemented using 2-D array for dense case and row-sparse for sparse case. Details are in a matrix series post.

Let m be a matrix, we could also use m.[i,j] to access the element, which is actually a function call to the the internal array. Thus there is some overhead. Some profiling is reported in a previous post.

The way to get rid of this kind of overhead is to program on the internal arrays directly. As in this answer in Stackoverflow.

Why F# is the language for data mining

2010-05-04T22:51:00.001+08:00

F# is a newly released language in Visual Studio 2010, as the 4-th language in .Net platform officially supported by Microsoft, aside with C#, VB.Net and C++/CLI. It has been a research project within Microsoft Research Cambridge for more than 5 years. Its origin could date back to the ML language designed by Turning Award winner Robin Milner in 1970s.

As a new yet old static-typed functional language designed for .Net platform, it will spread quickly in the coming years. It is such enjoyment to explore its usage in data mining, or more broadly, statistical computing area that I now share some my personal experience using it by trying to answer one question:

Why F# is the language for data mining.

This article is divided into 3 sections:

1. a review of data mining environments and a feature summary from them.

2. the answer. I list 7 reasons with 2 short examples.

3. a preview of a data mining library in F#, which will be open sourced this summer.

Before we go to the answer of this question, let’s see some existing data mining environments.

Data mining environments

Commercial:

SAS and SPSS are the main commercial data mining software. SAS has a mature macro language, also a reasonable GUI support. SPSS has a good GUI support. Their data mining algorithms are written in C/C++. The user doesn’t need to know about the underlying languages (C/C++), the data mining API is provided by a macro language (SAS) or GUI (SPSS). Using a macro language or GUI greatly eases users to perform data mining tasks.

There are also other products, e.g. the decision tree family developed by Salford.

Matlab is perhaps the most popular software in data mining community, especially in machine learning community. Although Matlab is not free, it is easily available in schools and research institutions. There are a lot of free programs in Matlab to perform data mining tasks. Matlab is so popular that putting an data mining algorithm name + "matlab" into Google gets what you need. The other reason is that the matrix operations in Matlab are efficient and convenient considering that the data tables in data mining are either dense or sparse matrices.

Free/Open source:

There are quite a lot of open source software of different quality. Here’s a list I am interested:

1. R. R is for statistics. The core part of R contains the R languages and the stable statistical procedures. The data mining part is a set of external packages. Most of the data mining models are available, e.g. Decision Trees, Neural Network, SVM, etc. R is a very capable language, however, it is not popular in data mining community. For data mining and machine learning, Matlab is more popular than R. R core does not support sparse data, thus a lot of extensions do not consider this support, which is extremely important. As in data mining, a lot of datasets are sparse, a tool that has sparse and large scale design in the beginning is very critical for performing data mining.

2. Weka. Weka is the most famous open source software in data mining and machine learning field. If data mining has three perspectives: database, machine learning and statistics[2]. Weka is from the second. Weka is widely used in data mining courses, in which instructors are able to show a lot of data mining algorithms to the students, turning the formulas in the book into several mouse clicks and the resulting figures and graphs. Who use Weka for serious data analysis? As I know, most people only use Weka’s decision tree implementation J48, because decision tree is really useful sometimes and Weka is handy to run to get the tree and its results are close to the standard C4.5's. But its usage is often on dense data set, as building a tree on sparse data sets costs a huge amount of memory in Weka. Weka does support sparse dataset, however, the performance is very bad, the speed is slow, the memory usage is unacceptable. The design in Weka is too general, thus sacrifices the speed.

Weka is written in Java. So the extensions to Weka, e.g. a different boosting scheme for trees, are usually written in Java. This kind of task is better done in a language like R or Matlab, i.e. a interpreted script language, not in Java. There are interpreted languages in JVM, most notably, Scala and Clojure. Calling Weka from Scala or Clojure will be a good idea.

3. Orange. Orange is written in Python and C++. All data mining algorithms are coded in C++ for performance’s sake. Python is used as the script language to call the algorithms and the language to build the GUI. Using Python in data mining is a great idea. It is interpreted, an concise language to read and write, easy binding to C/C++. Orange only supports dense data set.

4. Lush. Designed by Prof. Yann LeCun. He is a big figure in machine learning, a definite expert. There are not so many researchers, who write great programs, and write great papers at the same time. Prof. LeCun is one of the few. Lush is a Lisp environment with bindings to algorithms written in C/C++. The neural network in it is mature. Lisp is good, very expressive. However, nowadays, not many people write it.

after we have seen a set of data mining software, we can draw a conclusion on

the most favored features of a data mining environment:

1. Data mining algorithms are built in. The implementation should reliable and efficient. It should supports dense and sparse datasets equally well. It should also be robust, e.g. effectively handle missing values, large datasets, etc.

2. script language and iterative shell. The script language could be used to do data preprocessing, calling data mining procedures, visualization and report easily.

3. Math: Matrix and vector (both sparse and dense), linear algebra, common numerical optimizations. Sometimes, we need to design new algorithms or modified existing algorithms to suite one specific data mining application. In this case, these math tools are really important, e.g. a simple logistic regression could be done in several lines of code if these math modules are available.

4. GUI is preferred, but not a must. GUI is good for data exploration and visualization, e.g. Excel is used a front end for a lot of Business Intelligence solutions for its great GUI and stability. When performing data mining, most of the time, we don't need GUI (just as when we use Matlab, we seldom use its START menu and prefer command line).

5. A integrated system. You use tool A do download the data, use tool B to do clean the data, use tool C to do some matrix operations on the data, finally you get a dataset and put it into an SVM command line tool. Not the end, you still need to parse the ouput of the SVM to get the result. And think about you need to repeat this process for several times. What's wrong here? Tools are scattered. The Unix philosophy "do one thing and do it well" does not apply here because the data structure from one thing to the other is very complex. This is why we need a integrated system, at least data processing and model building are connected seamlessly, like in Matlab.

Weka has all the 5 points. However, none of them are very good. E.g. Java is hard to be used as a script language as you need to write a lot of lines of code and compile them to perform an easy task, and using Java also requires you to be very familiar with Weka's java classes, which is just huge. Java could also be used to do data preprocessing, however, few use it for this purpose. Python, Perl or F# are more appropriate.

F#

Following lists the notable features, which I think are helpful for data mining.

1. Succinct, static typing at the same time. These two features are combined together. When we write programs, we need to think about the programming logic and the programs, e.g. write the loop, define a variable, write the type for it, etc. Programming logic is the core thing, but every language has a syntax to obey, so we need to write a lot of other things, which could be distractive during programming because you need to think about more things at the same time. More succinct means you don’t need to write so many `declarations’ for your programs. Python is also succinct, however it is not static typing, which means typing checking is done during runtime. Dynamic typing is slow. In a lot of cases, F# does not need to write types, its other syntax is also succinct, it is also static typing at the same time. Combining these two together means we can focus on the programming logic and at the same time keep the performance.

2. Functional and programming style. F# is functional, thus the A, B, C of a functional programming also applies to F#. Here A could be ‘function as the first class or high order functions’, which means you could write:

a |> split |> Seq.filter notStopWord |> Seq.map stemmer

to perform string split, filtering using a self defined filter and doing word stemming in a single line!

In an OOP language, this would require to use a pipe line design pattern, which is far more complex than high order functions here.

3. Functional and parallel. A function is pure if it does not change anything outside the function, and every time you pass the same parameter into it, it gets the same result (this is a weaker version of purity). Based on this, we can easily parallel pure functions. In data mining, most of cases, we can simply split the data into several parts to get parallel execution.

As an example, following is a piece of code to perform multiple k-means clusterings at the same time:

if isParallel then
let dorepeats =
{0..restart-1} |> PSeq.map (fun _ -> doClustering(ds, k, maxiter, init, algo) :?> kmeansResult)
dorepeats |> PSeq.toList |> List.min :> IClusterResult
(* use async
let dore =
{0..restart-1} |> Seq.map (fun _ -> async { return doClustering(ds, k, maxiter, init, algo) :?> kmeansResult })
dore |> Async.Parallel |> Async.RunSynchronously |> Array.min :> IClusterResult
*)
else
//printfn "no parallel"
let dorepeats =
{0..restart-1} |> Seq.map (fun _ -> doClustering(ds, k, maxiter, init, algo) :?> kmeansResult)
dorepeats |> Seq.min :> IClusterResult

Note that the difference between sequential code and parallel code is minor.

4. OOP. At a high level of designing an application, the design patterns come from the application itself. A data mining application has its own working patterns. Design patterns are what OOP is good at. F# provides limited, but concise support for OOP. E.g. the following is the interface for a learner and predictor:

type IPredictor =
abstract member Predict: gvector -> int
abstract member Predict: datavecs -> int array
abstract member Predict: dataset2 -> float * (int array)

type IProbPredictor =
abstract member Predict: gvector -> int * (float array)
abstract member Predict: datavecs -> (int array) * (float array array)
abstract member Predict: dataset2 -> float * (int array) * (float array array)

type ILearner =
abstract member Learn: dataset2 * parameters -> IPredictor

A classification algorithm has two components, a learner and predictor. All learners implement ILearner interface. All predictors implement the IPredictor interface,some predictors also implement the IProbPredictor interface if it supports probability output. Because the code is so concise, I can change this interface more easily than in C#.

5. Interpret language. This is simply a must! Everyone wants to write a single line to perform data loading and model building together and see the result immediately.

6. .Net! If F# has the above 5, but without this one, it will not be successful. Tied with .Net has a lot of advantages, most notably a) a stable environment and a lot of stable libraries, b) writing components in C# when needed, c) the components in F# could be used by other .Net languages with no effort d) Visual Studio 2008 and 2010 are very good IDEs and finally e) Microsoft supports it, the `evil company’ is more responsible for its products than others and at last f) the performance of .Net is very good.

7. IO and preprocessing! This point is subsumed in previous points. However it is so important that I separate it out. A lot of software simply assume that data mining is the processing applying an algorithm to a dataset. This is WRONG! It is said that 80% of the data mining work is to preparing the data to get a dataset[1]. To prepare the data, we need to do a lot of IOs, from Internet, databases, xml files etc. .Net provides the most stable yet simple way to do these tasks. F# language is also well suited to perform data preprocessing. As an interpreted language, we can test and see the result immediately.

I use two examples from the code in my data mining library as examples to show F#'s strength.

Example 1. Wrod Tokenizer.

The following code is the Tokenizer module for text mining. It is so short, yet can perform various tokenizing tasks, removing stop words, doing stemming, bi-gram etc.. The conciseness comes from F# as a functional language and .Net library.

module Tokenizer =
let sep = [|' '; '\r'; '\n'; '-'; '.'; ',' ; '\t'; '!'; '?'; '\'' |]

let unigramToken (s:string) =
s.Split(sep, StringSplitOptions.RemoveEmptyEntries) |> Seq.filter (fun w -> not (isStopWord w)) |> Seq.map (fun s -> s.ToLower())

let unigramTokenWithStem s = s |> unigramToken |> Seq.map stemmer

let bigramToken (s:string) =
let uni = s.Split(sep, StringSplitOptions.RemoveEmptyEntries)
let bi = seq { for i=0 to uni.Length-2 do yield (uni.[i]+uni.[i+1]) }
bi |> Seq.filter (fun w -> not (isStopWord w)) |> Seq.map (fun s -> s.ToLower())

let bigramTokenWithStem s = s |> bigramToken |> Seq.map stemmer

let defaultToken = unigramTokenWithStem

Example 2. Parallel Tasks

You probably see similar code before as an example of the F# asynchronous programming. Here I changed the 'example' code to let it disguises as a Mozilla browser(this is useful when you try to crawling Google search results), also it does some retries and support exception handling. The code is so short!! Using it downloading 10000 posts (I need some special posts from a question answer forum for my research) in parallel costs only half an hour.

let fetchUrl (url:string) =
async {
try
let req = WebRequest.Create(url) :?> HttpWebRequest
req.UserAgent <- "Mozilla/4.0 (compatible; MSIE 6.0; Windows NT 5.1; SV1; .NET CLR 1.1.4322)"; req.Method <- "GET"; req.AllowAutoRedirect <- true;

req.MaximumAutomaticRedirections <- 4;

let! response1 = req.AsyncGetResponse() let response = response1 :?> HttpWebResponse
use stream = response.GetResponseStream()
use streamreader = new System.IO.StreamReader(stream)
return! streamreader.AsyncReadToEnd() // .ReadToEnd()
with
_ -> return "" // if there's any exception, just return an empty string
// checked and found that pages that cannot be retrieved mainly because the pages are already deleted
// not because of the downloader
}

let getPages pageIds rooturl =
pageIds |> Seq.map (fun id -> let name = rooturl + id.ToString() in fetchUrl name) |> Async.Parallel |> Async.RunSynchronously

A data mining library in F#

In the past month, I have written a prototype of an incomplete data mining system in F#. The more I write, the more I feel that it is the language for data mining. The library will be open sourced this summer after functionality completion and testing. I’d like to pronounce some highlights here:

1. Datasets and IO.

Dense and Sparse datasets are equally support! Various data formats are supported, csv files, weka's ARFF format, libsvm/svm light data format, Excel table, for both read and save. The library does not have its own data format. Processed data are stored in different formats for different needs. Using a built-in database will be considered in the future, e.g. the light-weight SQLite.

There are a set of preprocessing functions for datasets, e.g., filling missing values, feature value normalization, etc.

2. Data mining algorithms.

A few State-of-the-art data mining algorithms are implemented:

classification: decision tree family, logistic regression, support vector machines (by compiling libsvm to .Net using C++/CLI), neural networks, Naive Bayes, KNN.

regression: simple least square linear regression, neural networks and SVM(using libsvm).

clustering: k-means and its variants, EM on Gaussian Mixture Model, spectral clustering

dimension reduction: SOM, PCA, LDA.

feature selection: Mutual Information, decision trees, logistic regression with Wald test.

MISC: Probability Latent Semantic Analysis(PLSA), co-clustering

association rule and sequence mining: I have a plan in mind, but the actual implementation may delay.

The library is mainly for data mining, not for machine learning research. So only the most useful models are implemented. I hope after several years, the they will be very stable with few bugs. Trying to make it big would be dangerous. E.g. Weka is big in the sense of its algorithms, however, the standard statistical and data mining workhouse, logistic regression, is badly implemented there. If even logistic regression cannot apply to a modest text classification problem, what shall we expect Weka's other algorithms?

Also, implement a lot of immature models is less meaningful for practical data mining. E.g. the data mining module in SQL Server 2008 only implements a few state-of-the-art models. The novel algorithms are suitable to be implemented on demand for a specific application or dataset.

3. Parallel.

Some algorithms have parallel implementations in algorithm level, e.g. k-means clustering. At a higher level, I define a set of tasks using the type systems in F#, e.g. the classification task is defined as

type classificationEvaluationTask =
| CrossValidation of int * dataset2 * ILearner * parameters
| RandomSplit of float * dataset2 * ILearner * parameters
| TrainTest of dataset2 * dataset2 * ILearner * parameters
| LeaveOneOut of dataset2 * ILearner * parameters

Data mining tasks in the library run in parallel on default, e.g. performing a classification task and a clustering task at the same time, or, doing cross validation on different parameters to choose the best parameter of a model at the same time. A major technique to boost the performance of classifiers is ensemble learning, which bags a lot of models together. Parallel really helps a lot in this kind of learning task. Of course, parallel programming could be done in C++/Java using threads too. The advantage of using F# is that the interface is very easy to use, parallel executing various tasks at the same time needs only a single line.

4. A service provider design pattern for data mining algorithms. An algorithm, e.g. KNN will have several implementations, F# implementation is called on default, the user could also start a service provider and call the implementation of that provider. A candidate for such a wrapper is OpenCV’s machine learning component for its robustness, speed and friendly license.

ML library in OpenCV only supports dense dataset as machine learning tasks in CV usually use dense data format. I haven't found a good service provider for sparse data sets yet. But for large scale sparse data sets, the common practice is to use SVM, Logistic regression and tree family, which I already have code in C/C++.

5. Visualization. A limited support is already done by extending Luca Bolognese' code snippet. The visualization is mainly for understanding the data and diagnosing the algorithms, e.g. plotting the scatter of two columns of a data table, drawing the objective values through iterations. For high quality visualization, choose Mathematica, Matlab, R or Excel:)

6. Math library. It is a pity that there's no stable, yet free math library in .Net. For data mining, a math library is really important. I use the matrix and vector implementation in F#'s PowerPack. The matrix type in PowerPack supports both dense and sparse matrix, however sparse vector is not implemented. So I implemented the sparse vector myself. PowerPack now deprecates its math provider, which I pick up and use as the linear algebra module. So actually, I made little effort to get matrix supported.

Other math functions are implemented on demand. Some gradient based optimization solvers are implemented to solve linear regression and logistic regression, e.g. conjugate gradient decedent and LBFGS, which are two state-of-the-art gradient based solvers. Currently I implemented them in F#, the performance is modest compared to my C++ implementation. Code tuning is yet to be done. Some probability distributions are implemented for EM and Naive Bayes. Some important statistical algorithms, e.g., various tests, are also considered to be implemented in the future.

Acknowledgment
I'd like to thank Haipeng Zhang for reviewing this article for me.

[1] S. Zhang, C. Zhang, and Q. Yang. Data preparation for data mining. Applied Artificial Intelligence, 2003, 17(5-6): 375-381.

[2] Z.-H. Zhou. Three perspectives of data mining. Artificial Intelligence, 2003, 143(1): 139-146.

Matrix and linear algebra in F#, Part V: Sparse matrix implementation in PowerPack, and PInvoke a large scale SVD library as an application.

2010-04-10T14:34:00.001+08:00

In the previous posts of this series. We mainly talked about dense matrices in F#:

Part I: a tutorial on using matrix type in F#, mainly focused on dense matrix.

Part II: doing dense matrix linear algebra using math providers.

Part III: eigenface application

Part IV: performance issues

In this post, we move to sparse matrices. In many applications in data mining, the data table is stored as a sparse matrix. For example, in text mining, a document can be represented by a row-vector(or column-vector), the value at index i means that how many times word i occurs in this document. While the vocabulary could be quite large (tens of thousands), thus the document vector is highly sparse. If we store the documents in a matrix, then the matrix is also highly sparse. In DNA sequence analysis applications, the same situation occurs, that is we need to deal with sparse matrices!

We first dig into the implementation of F# sparse matrix in PowerPack. This is important, as we need to know how the sparse matrices are represented in F#, what algorithms are supported, etc. Based on the F# implementation, we can build our own routines, e.g. change the row based sparse representation to column based. At last, we PInvoke a SVD library in C as a example to show how to link existing stable linear algebra routines for sparse matrices.

The sparse matrix implementation in PowerPack

Sparse matrices have several different representation, which also have different strength. A good introduction is this Wikipedia page.

F# uses Yale Sparse Matrix Format, or another name, Compressed Sparse Row. (ref to the above Wiki article.)

Let’s check this out:

> let B = Matrix.initSparse 4 3 [0,0,2.3; 2,0,3.8; 1,1,1.3; 0,2,4.2; 1,2,2.2; 2,2,0.5]
val B : matrix = matrix [[2.3; 0.0; 4.2]
                     [0.0; 1.3; 2.2]
                     [3.8; 0.0; 0.5]
                     [0.0; 0.0; 0.0]]
> B.InternalSparseValues
val it : float [] = [|2.3; 4.2; 1.3; 2.2; 3.8; 0.5|]
> B.InternalSparseColumnValues
val it : int [] = [|0; 2; 1; 2; 0; 2|]
> B.InternalSparseRowOffsets
val it : int [] = [|0; 2; 4; 6; 6|]

Non-zero entries are sequentially in B.InternalSparseValues in a row-by-row fashion. B.InternalSparseColumnValues stores their column indices. B.InternalSparseRowOffsets stores the row offsets, for row i, it starts at B.InternalSparseRowOffsets[i] and ends at B.InternalSparseRowOffsets[i+1]-1 in the value array(B.InternalSparseValues).

This representation allows to take out one row of a sparse matrix very efficiently. In data mining, data instances/samples are usually stored as a row vector in a matrix, so it supports efficient handling on the instance level.

However, this representation does not support efficient column operations. In numerical computing, taking columns is more often (following the old Fortran tradition), thus a lot of libraries use Compressed Sparse Column representation for sparse matrices. It is also called Harwell-Boeing sparse matrix format (HB format for short) in the numerical computing community. An example is here.

PowerPack’s aim for providing the Internal* member functions of matrix type is to let users to interpolate with other libraries. However, it does not provide a function to convert the row representation to the column representation. Thus, this could be our first exercise:

let toColumnSparse (A:matrix) =
if A.IsDense then failwith "only for sparse matrix!"
let nrow = A.NumRows
let ncol = A.NumCols
let cidx = A.InternalSparseColumnValues
let roffset = A.InternalSparseRowOffsets
let vals = A.InternalSparseValues
let nval = vals.Length

// 1. scan to get column offset
let colCnt = Array.create ncol 0
for i=0 to nrow-1 do
    for idx=roffset.[i] to roffset.[i+1]-1 do
        colCnt.[cidx.[idx]] <- colCnt.[cidx.[idx]] + 1
let colOffset = Array.zeroCreate (ncol+1)
for i=1 to ncol do
    colOffset.[i] <- colOffset.[i-1] + colCnt.[i-1]
Array.Clear(colCnt, 0, ncol) // clear cnt array

// 2. score the value
let rowIdx = Array.create nval 0
let vals2 = Array.create nval 0.0
for i=0 to nrow-1 do
    for idx=roffset.[i] to roffset.[i+1]-1 do
        let j = cidx.[idx]
        vals2.[colOffset.[j] + colCnt.[j]] <- vals.[idx]
        rowIdx.[colOffset.[j] + colCnt.[j]] <- i
        colCnt.[j] <- colCnt.[j] + 1
    
vals2, rowIdx, colOffset

This implementation is quite efficient. It does two scans and returns three arrays for the HB format. This implementation currently only supports sparse matrices, it should also support dense matrices as we often need to transform a dense matrix into a sparse one. We can add the support in this function, or we can first transfer a dense matrix into a sparse one and use the above function. Both are OK. However, PowerPack currently does not support dense to sparse operation (it supports sparse to dense, which is used more often). Let this dense-to-sparse operation be our second exercise:

let toSparse (A:matrix) =

    if A.IsSparse then failwith "A should be a desne matrix"
let l =
    seq {
        let B = A.InternalDenseValues // use the internal array
        for i=0 to (Array2D.length1 B)-1 do
            for j=0 to (Array2D.length2 B)-1 do
                if B.[i,j] > 1e-30 then
                    yield (i,j,B.[i,j])
    }
Matrix.initSparse A.NumRows A.NumCols l

To get the HB format of a dense matrix A, we simply use:

let vals, rowIdx, colOffset = A |> toSparse |> toColumnSparse

At last of this section, let’s have a look at the annotated implementation for Matrix.initSparse in matrix.fs to have a better understanding of the F# sparse matrix type:

/// Create a matrix from a sparse sequence
let initSparseMatrixGU maxi maxj ops s =
(* 1.  the matrix is an array of dictionarys
       tab[i] is the dictionary for row i *)
let tab = Array.create maxi null
let count = ref 0
for (i,j,v) in s do
    if i < 0 || i >= maxi || j <0 || j >= maxj then failwith "initial value out of range";
    count := !count + 1;
    let tab2 =
        match tab.[i] with
        | null ->
            let tab2 = new Dictionary<_,_>(3)
            tab.[i] <- tab2;
            tab2
        | tab2 -> tab2
    tab2.[j] <- v
// 2. calcuate the offset for each row
// need to be optimized!
let offsA =
   let rowsAcc = Array.zeroCreate (maxi + 1)
   let mutable acc = 0
   // 3. this loop could be optimized using
   // sorted map for tab
   for i = 0 to maxi-1 do
      rowsAcc.[i] <- acc;
      acc <- match tab.[i] with
              | null -> acc
              | tab2 -> acc+tab2.Count
   rowsAcc.[maxi] <- acc;
   rowsAcc
// 4. get the column indices and values
let colsA,valsA =
   let colsAcc = new ResizeArray<_>(!count)
   let valsAcc = new ResizeArray<_>(!count)
   for i = 0 to maxi-1 do
      match tab.[i] with
      | null -> ()
      | tab2 -> tab2 |> Seq.toArray |> Array.sortBy (fun kvp -> kvp.Key) |> Array.iter (fun kvp -> colsAcc.Add(kvp.Key); valsAcc.Add(kvp.Value));
   colsAcc.ToArray(), valsAcc.ToArray()
// 5. call the SparseMatrix constructor
SparseMatrix(opsData=ops, sparseValues=valsA, sparseRowOffsets=offsA, ncols=maxj, columnValues=colsA)

As noted in the code, one possible optimization is to use SortedMap for table, rather than an array. But this SortedMap has more overhead, the current implementation is already good. The other possible way is to sort the (i,j,val) sequence, which avoids the overhead in using a Dictionary structure.

Only a few matrix operations are implemented for sparse matrices, e.g. +, – and * are supported. However, map, columns and rows are not supported. This does not quite matter as when we need sparse matrices, we will be usually dealing with large datasets. For large datasets, calling a specialized library or writing the code ourselves is a better solution, as we will see the SVD example blow:

PInvoke SVDLIBC

In a previous post, We already know how to write a simple matrix multiplication in C, and call it from F# using P/Invoke. Here we move to a more useful one, a large scale SVD library, SVDLIBC. For a small dense SVD, using lapack’s svd is just fine. However, for a 10000-by-10000 sparse matrix, we need a more powerful one. (ARPACK project is dedicated to this kind of decompositions. SVDLIBC is a C translation of a small part ARPACK.)

The SVDLIBC is a very good svd solver. It also provides a command line tool to do SVD for sparse or dense matrices. However, it uses some non-standard headers for I/O. To make it compile, we need to delete some code for IO. The main svd solver (in las2.c) is:

SVDRec svdLAS2A(SMat A, longdimensions)

A wrapper with a clear interface is need for this solver:

#define CHECK(ptr)  if (!(ptr)) return 0;

__declspec(dllexport)
int svds(int nrow, int ncol, int nval, double *val, int *row, int *offset,
     int dim,
     double *Uval,
     double *Sval,
     double *Vval)
{
SMat A;
SVDRec res;
DMat U, V; double *S;
FILE *fp;


A = svdNewSMat(nrow, ncol, nval);
CHECK(A);
memcpy(A->value, val, sizeof(double) * nval);
memcpy(A->rowind, row, sizeof(int) * nval);
memcpy(A->pointr, offset, sizeof(int) * (ncol+1));


res = svdLAS2A(A, dim);
CHECK(res);
CHECK(res->d == dim); // the dimension passed in must be correct!

S = res->S;
memcpy(Sval, S, sizeof(double) * dim);

U = svdTransposeD(res->Ut);    CHECK(U);
V = svdTransposeD(res->Vt);    CHECK(V);
memcpy(Uval, &U->value[0][0], sizeof(double) * (U->rows * U->cols));
memcpy(Vval, &V->value[0][0], sizeof(double) * (V->rows * U->cols));
svdFreeDMat(U);
svdFreeDMat(V);
svdFreeSMat(A);
svdFreeSVDRec(res);
return 1; // successful
}
and in F#:

module Native =
[<System.Runtime.InteropServices.DllImport(@"svdlibc.dll",EntryPoint="svds")>]
extern int svds(int nrow, int ncol, int nval, double *vals, int *row, int *offset, int dim, double *Uval, double *Sval, double *Vval);

module LinearAlgebra =
// Sparse SVD
// A: F# sparse matrix
// di: the dimension
let svds (A:matrix) (di:int)=
    /// A =  U * S * Vt
    if A.IsDense then failwith "only for sparse matrix!"
    else
        let nrow = A.NumRows
        let ncol = A.NumCols
        let nval = A.InternalSparseValues.Length
        // let dim = min nrow ncol
        let dim = max 1 (min (min nrow ncol) di) // choose a valid value
        let U = Matrix.zero nrow dim
        let S = Vector.zero dim
        let V = Matrix.zero ncol dim
    
        let colVals, rowIdx, colOffset = MatrixUtility.toColumnSparse A
        let valsP = NativeUtilities.pinA colVals
        let rowP = NativeUtilities.pinA rowIdx
        let offsetP = NativeUtilities.pinA colOffset
        let Up, Vp = NativeUtilities.pinMM U V
        let Sp = NativeUtilities.pinV S
    
        let ret = Native.svds(nrow, ncol, nval, valsP.Ptr, rowP.Ptr, offsetP.Ptr, dim, Up.Ptr, Sp.Ptr, Vp.Ptr)
        valsP.Free()
        rowP.Free()
        offsetP.Free()
        Up.Free()
        Vp.Free()
        Sp.Free()
        if ret = 0 then
            failwith "error in pinvoke svds"
        U, S, V
// default Sparse SVD     
// A: F# sparse matrix
let svds0 A =
    svds A (min A.NumCols A.NumRows)

test it:

let test2() =
let r = new System.Random()
let nrow = 1000
let ncol = 1000
let nval = 100000

let l =
    seq {
        for i=1 to nval do
            yield (r.Next()%nrow, r.Next()%ncol, r.NextDouble())
    }
let A = Matrix.initSparse nrow ncol l
tic()
let U, S, V = svds A 100
toc("svds")
()

We will in later posts to show the data mining applications of SVD.

NativeUtilities module

I used utility functions in NativeUtilities to convert the F# array/matrix/vector and their native array representations. Here’s its implementation from F# math-provider source code:

open System
open System.Runtime.InteropServices
open Microsoft.FSharp.NativeInterop
open Microsoft.FSharp.Math


// from math-provider source code
module NativeUtilities = begin
let nativeArray_as_CMatrix_colvec (arr: 'T NativeArray) = new CMatrix<_>(arr.Ptr,arr.Length,1)
let nativeArray_as_FortranMatrix_colvec (arr: 'T NativeArray) = new FortranMatrix<_>(arr.Ptr,arr.Length,1)
let pinM m = PinnedArray2.of_matrix(m)
let pinV v = PinnedArray.of_vector(v)
let pinA arr = PinnedArray.of_array(arr)

let pinA2 arr = PinnedArray2.of_array2D(arr)

let pinMV m1 v2 = pinM m1,pinV v2
let pinVV v1 v2 = pinV v1,pinV v2
let pinAA v1 v2 = pinA v1,pinA v2
let pinMVV m1 v2 m3 = pinM m1,pinV v2,pinV m3
let pinMM m1 m2  = pinM m1,pinM m2
let pinMMM m1 m2 m3 = pinM m1,pinM m2,pinM m3
let freeM (pA: PinnedArray2<'T>) = pA.Free()
let freeV (pA: PinnedArray<'T>) = pA.Free()
let freeA (pA: PinnedArray<'T>) = pA.Free()

let freeA2 a = freeM a

let freeMV (pA: PinnedArray2<'T>,pB : PinnedArray<'T>) = pA.Free(); pB.Free()
let freeVV (pA: PinnedArray<'T>,pB : PinnedArray<'T>) = pA.Free(); pB.Free()
let freeAA (pA: PinnedArray<'T>,pB : PinnedArray<'T>) = pA.Free(); pB.Free()
let freeMM (pA: PinnedArray2<'T>,(pB: PinnedArray2<'T>)) = pA.Free();pB.Free()
let freeMMM (pA: PinnedArray2<'T>,(pB: PinnedArray2<'T>),(pC: PinnedArray2<'T>)) = pA.Free();pB.Free();pC.Free()
let freeMVV (pA: PinnedArray2<'T>,(pB: PinnedArray<'T>),(pC: PinnedArray<'T>)) = pA.Free();pB.Free();pC.Free()

let matrixDims (m:Matrix<_>) = m.NumRows, m.NumCols
let matrixDim1 (m:Matrix<_>) = m.NumRows
let matrixDim2 (m:Matrix<_>) = m.NumCols
let vectorDim  (v:Vector<_>) = v.Length

let assertDimensions functionName (aName,bName) (a,b) =
  if a=b then () else
  failwith (sprintf "Require %s = %s, but %s = %d and %s = %d in %s" aName bName aName a bName b functionName)
end

A Note on Compressed Sparse Row and Compressed Sparse Column(HB format)

We can see both representations have deficiency, is there a magic structure or an engraining trick taking the advantage of both? Let’s check the standard software Matlab.

In Matlab, we usually write A(:,j) to take the j-th column of sparse matrix A or A(i,:) to take the i-th row. Does Matlab have a magic to let both operations run efficiently. Then answer is NO. The following script is used to test this:

function [ ret ] = testOctaveSvds ()

nrow = 10000;
ncol = 10000;
nval = 100000;
rows = floor(rand(1, nval)*nrow)+1;
cols = floor(rand(1, nval)*ncol)+1;
vals = rand(1, nval);
m = sparse(rows, cols, vals, nrow, ncol);
tic;
s = 0;
for i=1:10000,
c = floor(rand(1) * ncol) + 1;
s = s + sum(m(:,c));
end
fprintf('cols: ')
toc;

tic;
s = 0;
for i=1:10000,
r = floor(rand(1) * nrow) + 1;
s = s + sum(m(r,:));
end
fprintf('rows: ')
toc;

endfunction

In Matlab 2009a, the output is:

cols: Elapsed time is 0.204693 seconds.

rows: Elapsed time is 16.058717 seconds.

Octave also has similar result. From the result, we can deduce that Matlab/Octave use HB format for sparse matrix and it does not do heavy optimization for the operation for taking the rows. This is reasonable, as mentioned before, that when using sparse matrices, the user/programmer has the obligation to optimize the program, rather than the matrix library.

PInvoke native C/C++ libraries in F#

2010-04-05T01:00:00.001+08:00

In this blog post, we explore on PInvoke in F#. Before going into PInvoke, an matrix multiplication example is also given to show how to build a DLL for PInvoke, which is important as in most of the cases we don't have pre-built DLLs.

PInvoke stands for Platform Invoke, meaning calling naive libraries. The two major reasons to use PInvoke is 1) native performance, native code is usually faster than .Net managed code, especially for numerical computing routines, and 2) existing stable(well tested) libraries can be used on the fly on .Net.

Of course, it also has some disadvantages, e.g. 1) extra instructions are used for each PInvoke call, this is not important as long as P/Invoke is not used too frequently, 2) calling native code means not safe anymore, e.g. no exception, no memory access checking.

There are several tutorials online for C#: here, here and here. However, these are more focused on existing windows dlls, their motivation for PInvoke is that some functions are not available on .Net, thus calling an existing native DLL, e.g. Win32 API. However, we are talking about performance in this series, PInvoke means that we rewrite the performance critical part in C/C++, or compile an existing library (maybe we also need to make a C wrapper of it in the meantime) and call it using PInvoke. So we should start with the native code part as we don’t have the native dll yet.

Matrix multiplication

Because C and C++ are slightly different, and C++ usually compiles C code well. So we’ll deal with C++ mode in VC++ compiler, i.e. all files have .cpp extensions.

First we write mattrixproduct.cpp:

__declspec(dllexport) int minus(int a, int b)
{
return a - b;
}

extern "C" __declspec(dllexport) int add(int a, int b)
{
return a + b;
}

extern "C" __declspec(dllexport)
void matmul(double *a, int an, int am,
        double *b, int bn, int bm,
        double *c)
{
int i, j, k;

for (i=0; i<an; i++)
    for (j=0; j<bm; j++) {
        double s = 0;
        for (k=0; k<am; k++)
            s += a[i*am + k] * b[k*bm + j];
        c[i*bm+j] = s;
    }
}

__declspec(dllexport) is the header we must add to each function to tell the compiler to generate this function in the DLL. extern "C" is to tell the compiler that this is a C code. We use command line:
cl.exe /LD mattrixproduct.cpp
to get mattrixproduct.dll. We can also create a new Visual C++ project and set the target to .dll.

Here’s the information for this dll: (use command line: link /DUMP /EXPORTS mattrixproduct.dll)

ordinal hint RVA name
1 0 00001000 ?minus@@YAHHH@Z = ?minus@@YAHHH@Z (int __cdecl minus(int,int))
2 1 00001010 add = _add
3 2 00001020 matmul = _matmul

We can see that minus function (without extern “C”) has a non-readable name, while the other two have correct names.

Now, we have the DLL, let’s move on to the PInvoke in F#. Here is the code:

open System
open System.Runtime.InteropServices
open Microsoft.FSharp.NativeInterop
open Microsoft.FSharp.Math

module Native =
[<System.Runtime.InteropServices.DllImport(@"mattrixproduct.dll",
EntryPoint="add")>]
extern int add(int a, int b);

// notice the wired name for minus
[<System.Runtime.InteropServices.DllImport(@"mattrixproduct.dll",
EntryPoint="?minus@@YAHHH@Z")>]
extern int minus(int a, int b);

[<System.Runtime.InteropServices.DllImport(@"mattrixproduct.dll",
EntryPoint="matmul")>]
extern void matmul(double *a, int an, int am, double *b, int bn,
int bm, double *c);


let a = Native.add(10,20)
let b = Native.minus(10,20)
let A = matrix [[1.0;2.0;]; [3.0;4.0;]]
let B = matrix [[3.0;2.0;1.0]; [1.0;1.0;1.0]]
let C = Matrix.zero 2 3
let Ap = PinnedArray2.of_matrix(A)
let Bp = PinnedArray2.of_matrix(B)
let Cp = PinnedArray2.of_matrix(C)
Native.matmul(Ap.Ptr, 2, 2, Bp.Ptr, 2, 3, Cp.Ptr)
Ap.Free()
Bp.Free()
Cp.Free()

printfn "a = %A\nb = %A\nC = %A" a b C

There are several things to mention:

1. For basic type, like int, double, as in the add and minus, they don’t need any special treatment.

2. The hard part is the array. This introduces concept called data marshaling. Because .Net arrays and native arrays are different, so they cannot be directly mapped. This is why we open namespaces for Interop. Another thing is that data marshaling does not support 2 dimensional arrays directly! This is why I used 1-D array in the C implementation. High performance code usually only use 1-D arrays, so this inconvenience does not quite matter to us. To marshal the array, we get a one dimensional pointer for a matrix using PinnedArray2.of_matrix, which is defined in F# PowerPack (in NativeArrayExtensions.fs), and then pass this pointer in the function call.

The performance

Let's now test the performance of the wrapped native matrix multiplication. First write some test on both small matrices and big matrices:

let mm A B =
let Ap = PinnedArray2.of_matrix(A)
let Bp = PinnedArray2.of_matrix(B)
let C = Matrix.zero (A.NumRows) (B.NumCols)
let Cp = PinnedArray2.of_matrix(C)
Native.matmul(Ap.Ptr, A.NumRows, A.NumCols, Bp.Ptr, B.NumRows, B.NumCols, Cp.Ptr)
Ap.Free()
Bp.Free()
Cp.Free()
C

let test() =
// use small matrices
let A = matrix [[1.0;2.0;]; [3.0;4.0;]]
let B = matrix [[3.0;2.0;1.0]; [1.0;1.0;1.0]]
tic()
mm A B |> ignore
toc("my native matrix *")
A * B  |> ignore
toc("F# *")

// use big matrices
let r = new System.Random()
let A1 = Matrix.init (92*112) 280 (fun i j -> r.NextDouble())
let A' = A1.Transpose
tic()
let C1 = mm A' A1
toc("my native matrix *")
let C2 = A' * A1
toc("F# *")

if C1.Equals(C2) then true
else false

The timing functions tic() and toc() are defined in Part IV of the matrix and linear algebra series. Run it in F# interactive:

> test();;
my native matrix * 1.758400 ms
F# * 0.005100 ms
my native matrix * 7207.073800 ms
F# * 16524.670200 ms
val it : bool = true

First the return value is true, indicating that my C implementation is correct. In the small matrix case, the overhead of PInvoke is obvious. In the big matrix case, the performance boost using native code is also obvious.

Note 1: C++ class

So far so good. Because we didn’t touch C++ yet! Classes, templates are more complex.

To work with these. One method is to write C wrappers for C++ classes. This method is commonly used in practice.

We can also directly add __declspec(dllexport) to the class definition. However, creating objects via constructors is not supported anymore. One pattern is to write two static methods: one for creating an object and the other for freeing an object. Member functions are used by providing this pointer explicitly. If you know how to enable object style programming in C, this patter would be familiar to you. Anyway, in this style, the accessing ability is equivalent to C’s. As this pattern requires modifying existing class definitions(although only a little), it is less commonly used, thus I don’t give detailed example here.

Note 2: Pure C

The above example is given in C++ as in most of the cases we are dealing with C++ files, at least C files could be compiled in C++.

For Pure C files (with .c extensions), things are actually easier. Just removing extern "C" would be OK.

Using existing C/C++ code without PInvoke

We can use C++/CLI to compile a C/C++ library without modifying anything into a unsafe managed dll. I will write a separate blog for this later.

Remark

In this tutorial, we know the basics of P/Invoke in F#. We don’t touch two things:

1) Marshaling complex parameters, e.g. different format of strings, non-array pointers, C++ classes, etc.

2) MEMORY! MEMORY! MEMORY! The hard and dangerous part of P/Invoke is memory management. I don’t have enough experience of this part yet. Currently I just follow the style in Math-Provider.

Matrix and linear algebra in F#, Part IV: profile your program, find the bottleneck and speed it up: using matrix multiplication as an example

2010-04-03T21:04:00.001+08:00

We will use the face recognition program we developed in Part III to illustrate how to profile programs.

The performance problem of the face recognition program is that it takes about 25 seconds to process 280 faces (each face has 10304=112*92 pixels). I want to find which part takes most of the time. Is it due to the eigen decomposition?

Profile

Profiling means knowing how much time and memory every part of your program take. Good profiling tools could give a detail report on the running behavior of a program. In .Net, there are several good tools for this purpose. Some of them are free.

However, these tools are too big to use in a explorative/interactive programming style. Following Matlab’s tic and toc, I have written them for F#:

let stopWatch = new System.Diagnostics.Stopwatch()

let tic() =
stopWatch.Reset()
stopWatch.Start()

let toc(msg:string) =
stopWatch.Stop()
printfn "%s %f ms" msg stopWatch.Elapsed.TotalMilliseconds
stopWatch.Reset()
stopWatch.Start()

I also use the same style for my C programs:

#include <time.h>
time_t stopWatch;
void tic()
{
stopWatch = clock();
}
void toc(const char *msg)
{
double s = (clock() - stopWatch) * 1000.0 / CLOCKS_PER_SEC;
printf("%s: %.2lf ms\n", msg, s);
stopWatch = clock();
}

This implementation enables us to use a sequence of tocs to profile a code block line by line without the need to write tic every time. However this simple profiler can not be nested as there is only a global timer.

After putting several tocs into my eigCov:

tic()
let colm = colMean faces
toc("mean face")
let A = Matrix.init nrow ncol (fun i j -> faces.[i,j] - colm.[i])
toc("minus mean face")
let L = A.Transpose * A
toc("get L")
let val1, vec1 = eig L 
toc("lapack eig")
let v = val1 * (1.0 / float (ncol - 1))
let u = A * vec1.Transpose
toc("get u")
// normalize eigen vectors
let mutable cnt = 0
for i=0 to ncol-1 do
if abs v.[i] < 1e-8 then
   u.[0..nrow-1,i..i] <- Matrix.create nrow 1 0.0
else
   cnt <- cnt + 1
   let norm = Vector.norm (u.[0..nrow-1,i..i].ToVector())
   u.[0..nrow-1,i..i] <- u.[0..nrow-1,i..i] * ( 1.0 / float norm )
toc("normolize u")

It gets the following output:

mean face 639.891700 ms
minus mean face 195.484600 ms
get L 12885.751500 ms
lapack eig 247.015700 ms
get u 7169.749100 ms
normolize u 623.919700 ms

The bottle neck is in "Get L” and “get u”, both of which are matrix multiplication operations. The main computation part “lapack eig” actually costs very little time (as the 280-by-280 matrix is not big).

The matrix multiplication in “get L” is 280x10304 multiplies 10304x280, which is quite large!

Matrix Multiplication

Matrix multiplication is a great example to show that why we need well tuned libraries to perform numerical computations.

Table: the time cost for self-made(not optimized) implementations

implementations	operator * for matrix type	3 F# loops(use matrix)	3 F# loops(use float Array2D)	3 C loops (native)
time(seconds)	13.07	38.81	13.52	8.10

I used the following 3 loops:

for (i=0; i<N; i++)
for (j=0; j<N; j++) {
   double s = 0;
   for (k=0; k<M; k++)
       s += a[i][k] * b[k][j];
   c[i][j] = s;
}

The two F# programs are similar to the above one in C. We can find that if we use a lot of F# matrix element access operator [,], the performance is very poor. The * operator of matrix type does not have any optimization in it, so it costs the same amount of time as 3-loops (Array2D) does. Native code is faster than the managed ones as it does not do boundary checking for index.

Matlab, R and Numpy/Python all wrap optimized matrix implementations. Let us next see how they perform:

Table: the time cost for optimized implementations

Software	Matlab 2009a	R 2.10.1	Numpy 1.3.0
time	0.25	1.98	0.35

Matlab is the fastest, but it uses 200% CPU. While R and Numpy are single threaded. R takes 2 seconds. Matlab uses Intel MKL, Numpy uses Lapack(maybe optimized version), I don’t know what algebra routines R uses.

We can see the difference between optimized versions and non-optimized ones is very large. A similar report from a .Net math software vender is here. The conclusion is that we should use optimized matrix multiplication procedure when matrices are large.

Add optimized matrix multiplication support to math-provider

The math-provider has made an wrapper for dgemm_, which is a BLAS procedure for computing:

C  :=
alpha*op( A )*op( B ) + beta*C
where op (A) = A or A’.

This is general case of matrix multiplication. The math-provider also has an wrapper for dgemv_, matrix-vector multiplication.

However, they are not directly callable from Lapack service provider.

add the following code to linear_algebra_service.fs:
let MM a b =

  Service().dgemm_(a,b)

and add the following to linear_algebra.fs:

let MM A B =
if HaveService() then
  LinearAlgebraService.MM A B
else
  A * B

and add the following to linear_algebra.fsi:
/// BLAS matrix multiplication

val MM: A:matrix -> B:matrix -> matrix

Now we define a local name for it:

let mm = LinearAlgebra.MM

Using BLAS matrix multiplication routine costs about 4 seconds, which is faster than native C’s performance(8 seconds), but worse than Matlab’s or Numpy’s. The reason is that the BLAS(Netlib BLAS 3.1.1) is not optimized for my platform. If I use ATLAS for BLAS and use an Intel Fortran Compiler, the performance will be close to Numpy’s or Matlab’s.

By using BLAS’s matrix multiplication routine for the face recognition, the total running reduces from 25 seconds to 9 seconds.

Attachment

wait.

Matrix and linear algebra in F#, Part III: Eigen decomposition and face recognition

2010-04-03T17:55:00.001+08:00

The theory

Face recognition is one of the most import research problems in computer vision. It is also an important application for everyday use. A lot of security system has a face recognition component, as well as other parts, e.g., finger print recognition.

While this application has a lot of algorithms, the most famous one must be Eigen face proposed in

M. Turk and A. Pentland (1991). "Eigenfaces for recognition". Journal of Cognitive Neuroscience 3 (1): 71–86.

The theory of Eigen face is very simple if one has learned the concept of Principal Component Analysis (PCA). Given a set of faces, each face is a n-by-m matrix, if we reshape the matrix into a vector, then each face can be viewed as a point in the (n*m) dimensional space. The idea of Eigen face is to find some directions, which are orthogonal, in this space such that the points(faces) have a minimal representation, or have maximum variance along these directions. As you can see, this is an optimization problem. The magic part is that this optimization is solved by Singular Value Decomposition (SVD).

The following is a summary of the procedures to get the eigen faces:

1. for a set of faces, F = {f1,…,fn}, calculate its mean face m. F is a d-by-n matrix, where d is the number of pixels of an face and n is the number of faces in the dataset. fi is a column vector.

2. subtract the mean face from F, get A = {f1-m, … fn-m}. We want to get the eigenvalues and eigenvectors of Cov(A) = A * A’. The eigenvectors associated with big eigenvalues are eigenfaces. (In the program part, we will visualize these eigenfaces).

However, Cov(A) is usually too big. (If there are 10K pixels in each face, then the size of Cov(A) would be 10K-by-10K). A trick is to calculate the eigenvalues and eigenvectors of A’ * A (which is of small size if the number of faces n is small.) and calculate Cov(A)’s based on them. This trick is well documented in the program section.

Once we get the Eigen faces, we can do face recognition. As said before, each eigen face is actually an direction, we can map our new face and the face in the database into the space defined by the set of Eigen faces. For example, if we have 10 Eigen faces, and we know that they are orthogonal to each other, thus a 10-D real space is defined by the 10 Eigen faces. Once real faces are mapped into this 10-D space, we can measure the distance between the new face and the existing faces, and find the closet one as the label face for the new face.

Here are some reference for interested readers:

1. Eigenface wiki page. It has the mathematics and links to programs and papers.

2. A Matab Tutorial. and this one. These two are very good tutorials using Matlab.

3. Principal Component Analysis (PCA), Singular Value Decomposition and Eigendecompsition. If you want to know more theory of it, you can move on to have a good understanding of PCA.

The ORL face dataset

There are more than 10 face datasets available online, which range from different complexity. I use an quite old one the ORL face dataset. This one is well pre-processed, e.g. under the same light condition, of the same size, etc.

Using this dataset greatly eases our work in processing the images and enables us to focus on the Eigen face algorithm in F#.

This dataset contains 400 faces from 40 subjects, with each having 10 faces. The images are of .PGM format. For the I/O of the images, I wrote a separate blog post on reading, writing and visualizing this image format. That post also features on how to write a simple GUI in F#.

The data structures for the dataset:

type pgm = matrix
type person = matrix array
type dataset = person array

Each PGM image is an F# matrix, each person is an array of matrices, and the dataset is an array of persons.

Then we define functions to read the dataset for a given folder:

let readPerson folder =
  seq {
      for i=1 to 10 do  
          let name = folder + string(i) + ".pgm"
          yield readPgm name
  } |> Seq.toArray

let readDataset folder =
  seq {
      for i=1 to 40 do
          let name = folder + "s" + string(i) + "\\"
          yield readPerson name
  } |> Seq.toArray

We also define functions to split the dataset into training data and testing data. Training data are used to get the eigenfaces, and then faces in the testing data are recognized using the eigenfaces. For each person, the first 7 faces are used for training, and the 3 left ones are used for testing.

let splitDataSet (ds:dataset) =
  let train = seq {
      for p in ds do
          yield (
              seq {
              for i=0 to 9 do
                  if i<7 then yield p.[i]
              } |> Seq.toArray )
      }
    
  let test = seq {
      for p in ds do
          yield (
              seq {
              for i=0 to 9 do
                  if i>=7 then yield p.[i]
              } |> Seq.toArray )
      }
  train |> Seq.toArray, test |> Seq.toArray

let strip (data: 'a array array) =
  let dat = seq {
      for x in data do
          for y in x do
              yield y
      }
  let label = seq {
      let i = ref 0
      for x in data do
          for y in x do
              yield !i
          i := !i + 1
      }
  dat |> Seq.toArray, label |> Seq.toArray

Function strip is used to transform a data set (matrix array array) into a flat array (matrix array) and remember their labels.

The program

First we write some help functions:

let mat2vec (A:matrix) =
  let n = A.NumRows
  let m = A.NumCols
  Vector.init (n*m) (fun i->A.[i%n, i/n])

let vec2mat (vec:vector) nrow ncol =
  if vec.Length <> nrow * ncol then failwith "vec length != nrow * ncol "
  Matrix.init nrow ncol (fun i j -> vec.[j*nrow+i])

let colSum (m:matrix) =
  let nrow = m.NumRows
  Matrix.foldByRow (fun a b -> a+b) (Vector.create nrow 0.0) m

let colMean (m:matrix) =
  let ncol = m.NumCols
  colSum m * (1.0 / float ncol)

mat2vec is concate the columns of a matrix into a long column vector, while vec2mat does the inverse. Their behaviors are similar to the reshape function in Matlab. colSum and colMean are similar to sum and mean in Matlab.

To use the eigen decomposition for a symmetric matrix, we call lapack via math-provider. Read Part II of this F# matrix series to know to how use them.

open Microsoft.FSharp.Math

open Microsoft.FSharp.Math.Experimental
let isSucc = Experimental.LinearAlgebra.Lapack.Start()
printfn "service = %A" (if isSucc then "Netlib lapack" else "Managed code")
let eig = LinearAlgebra.EigenSpectrumWhenSymmetric

eig is a shorthand for EigenSpectrumWhenSymmetric. Notice that eig gets the eigenvectors in a matrix, each row of which is an eigenvector. This may be a bug in the math-provider wrapper, but since it only deals with square and symmetric matrices, this incontinence does not hurt too much.

The main function:

let eigCov (faces:matrix) numVecs =

    // each column is a face

  let nrow = faces.NumRows
  let ncol = faces.NumCols
  let colm = colMean faces
  let A = Matrix.init nrow ncol (fun i j -> faces.[i,j] - colm.[i])
  // A - each column is a difference face
  // Cov(A) = A * A', a big matrix
  // trick: A'*A*v(i) = lambda(i)*v(i)
  //      A*A'*A*u(i) = lambda(i)*A*u(i)
  //      (A*u(i)) is an eignvector of Cov(A)
  // L(A) = A' * A
  let L = A.Transpose * A
  let val1, vec1 = eig L // each row of vec1 is an eigenvector
  // get the eigen values and eigen vectors for Cov(A)
  let v = val1 * (1.0 / float (ncol - 1))
  let u = A * vec1.Transpose

  // normalize eigen vectors
  let mutable cnt = 0
  for i=0 to ncol-1 do
      if abs v.[i] < 1e-8 then
          u.[0..nrow-1,i..i] <- Matrix.create nrow 1 0.0
      else
          cnt <- cnt + 1
          let norm = Vector.norm (u.[0..nrow-1,i..i].ToVector())
          u.[0..nrow-1,i..i] <- u.[0..nrow-1,i..i] * ( 1.0 / float norm )

  // sort
  // let u' = u.PermuteColumns (fun j -> j ) //u.NumCols - j - 1)
  // the PermuteColumns function is buggy
  let u' = Matrix.init u.NumRows u.NumCols (fun i j -> u.[i, u.NumCols - j - 1])

  let v' = v.Permute (fun j -> v.NumRows - j - 1)

  let select = min cnt numVecs
  //printfn "%A" select
  if cnt < numVecs then
      printfn "warning: numvec is %d; only %d exist." numVecs cnt

  let v1 = v'.[0..(select-1)]
  let u1 = u'.Columns(0, select)
  let transformed = (u1.Transpose) * A
  v1, u1, colm, transformed

and:

let eigenface (p:matrix array) numVecs =
  let B = p |> Array.map mat2vec
  let nrow = B.[0].Length // nrow is big  == # of elements in a face
  let ncol = B.Length  // # of faces for a person
  // put a 2d face into a vector as a column in matrix A
  let A = Matrix.init nrow ncol (fun i j -> B.[j].[i])

  let v, u, mface, transformed = eigCov A numVecs
  (*
  (new PGMViewer(vec2mat mface 112 92, "mean face")).Show()
  for i=0 to 2 do
      let eface = scale (u.Column(i))
      let f = vec2mat eface 112 92
      (new PGMViewer(f, "eigenface " + i.ToString())).Show()
  *)

  v, u, mface, transformed

This function transform an image matrix into a vector and call the eigCov to get the eigenfaces. The commented lines in this function is used to output the mean face and the first 3 eigenfaces:

Use eigenfaces to do recognition for a new face:

let recognition (u:matrix) (mface:vector) (transformed:matrix) (face:vector) =

    let d = face - mface
  let t = u.Transpose * d

  let distance a b = Vector.norm (a - b)
  {0..transformed.NumCols-1}
  |> Seq.map (fun j -> distance t (transformed.Column(j)), j)
  |> Seq.min
  |> snd

To put the recognition into a learning and testing framework:

let learning data =
  let train, test = splitDataSet data
  let dat, lab = strip train
  let dat0, lab0 = strip test
  let v, d, mface, tran = eigenface dat 10

  let correctCnt =
      dat0
      |> Seq.mapi (fun i f ->
              let idx = recognition d mface tran (mat2vec f)
              let label = lab.[idx]
              if label = lab0.[i] then 1 else 0
          )
      |> Seq.sum
  float correctCnt / (float dat0.Length)

Using 10 most significant eigen faces from the 280(7 * 40) faces, we recognize the remaining 120(3 * 40) faces, and get 113 faces recognized to the correct person, thus a 94.17% accuracy, which is reasonable good considering that we have done any normalization on the faces, and the classifier is only 1 nearest neighbor.

Remark

In this post, we have seen that how F# could be used as an prototype tool like Matlab. In the matrix computation part, we need to write some utility functions ourselves, while in Matlab they are built in. However, the strength of F# comes from the language part, we could control the training and testing more freely than in Matlab. Moreover, once we have developed this face recognition component, it is just straightforward to use C# or VB.Net to call it.

Although we have a high accuracy in the final test, the result is only for a well prepared dataset. In real situations, things would be very hard. Maybe over 90% time is spent on tuning best parameters and preprocessing the face images.

Attachment

wait.

A PGM Image Viewer in F#, with PGM writer and reader

2010-04-03T17:49:00.001+08:00

This post is to serve the face recognition project.

In that project, we need a little tool to visualize/display Eigen faces. It would be very helpful during debugging/testing to display the Eigen faces and mean faces visually after they are calculated. We would also like to display an image to see if our IO part is correct or not.

The ORL face data set I use is of PGM format, so I need to load these images and display them.

First let’s get familiar with the PGM image format. PGM stands for Netpbm grayscale image. Its format is every easy:

P5 or P2
# 0 or more lines for comments
nrows ncolumns
max-value
image content

The format of image content part could be either text(P2) or binary(P5). More detail is here and here.

Before reading the file, let’s define the data structure for a PGM image. Because it is grey scale, thus it is actually an F# matrix:

type pgm = matrix

We would use two readers: a BinaryReader for binary image content and a StreamReader for text image content. (I think we could also use StreamReader for both, But I started using BinaryReader for binary image content and added the text image content later, so two readers are used.)

exception PgmFormatError of string
let readPgm path =

   let rec skipCommentLine br =
     let l = getline br
     if l.StartsWith("#") then
         skipCommentLine br
     else
         l
 let br = new BinaryReader(File.Open(path, FileMode.Open))
 let pgmType = getline br
 let sizeStr = skipCommentLine br
 let s = sizeStr.Split ' '
 let ncol = int (s.[0])
 let nrow = int (s.[1])
 printfn "nrow = %d ncol = %d" nrow ncol
 let maxStr = getline br

 if (pgmType.StartsWith("P5")) then
     // binary format
     let mutable b = Array.create (nrow * ncol) 0uy
     let nread = br.Read(b, 0, (nrow*ncol))
     let bb = b
     Matrix.init nrow ncol (fun i j -> float (bb.[i*ncol+j]))
 elif (pgmType.StartsWith("P2")) then
     // text format
     br.Close() // close the file and re-open
     use sr = new StreamReader(path)
     sr.ReadLine() |> ignore
     let mutable line = sr.ReadLine()
     while line.StartsWith("#") do
         line <- sr.ReadLine()
     sr.ReadLine() |> ignore

     let str =
         seq {
             while not sr.EndOfStream do
                 yield sr.ReadLine()
         }
         |> Seq.fold (fun a b -> a + b) ""
      
     let pic = Matrix.create nrow ncol 0.0
     let vals = str.Split([|' ';'\t';'\r';'\n'|], StringSplitOptions.RemoveEmptyEntries)
     printfn "%A" vals
     Matrix.init nrow ncol (fun i j -> float (vals.[i*ncol+j]))
 else
     // unknown format
     raise (PgmFormatError(path))

The getline function for BinaryReader:

let getline (br:BinaryReader) =

   let chSeq = seq {
     let ch = ref (br.ReadChar())
     while !ch <> '\n' do
         yield !ch
         ch :=  br.ReadChar()
     }
      
 let chArr = chSeq |> Seq.toArray
 let line = new String(chArr)

 line

Next we write the writer for a pgm image, I only implemented the writer for the text content case as this is mainly used for outputting some intermediate results, thus need to be human-readable.

let writePgm (path:string) (p:pgm) =
use fs = new FileStream(path, FileMode.Create)

   use sw = new StreamWriter(fs, Encoding.ASCII)
 sw.Write("P2\n{0} {1}\n255\n", p.NumCols, p.NumRows)

 for i=0 to (p.NumRows - 1) do
     for j=0 to (p.NumCols - 1) do
         sw.Write(" " + string (p.[i,j]))
     sw.WriteLine()

Finally, the viewer. F# actually is quite good at doing small GUI programs. First inherits from a form:

type PGMViewer(pic:matrix, name:string) as this =

   inherit Form(Width = 300, Height = 300)

then create a main menu, and menu items in it, a PictureBox to display bitmap images and add a SaveFileDialog:

let mainmenu = new MainMenu()

let mFile = new MenuItem()
let miFileSave = new MenuItem()
let miExit = new MenuItem()
let pb = new PictureBox()
let filesave = new SaveFileDialog()

Set necessary properties of these objects (put these code in a do block)

do
 // set the property of picture box
 let nrow = pic.NumRows
 let ncol = pic.NumCols
 let image = new Bitmap(ncol, nrow)
 printfn "%A %A" nrow ncol
 for i=0 to nrow-1 do
     for j=0 to ncol-1 do
         image.SetPixel(j, i, Color.FromArgb(int pic.[i,j], int pic.[i,j], int pic.[i,j]))
 pb.Image <- image
 pb.Height <- nrow
 pb.Width <- ncol

 // add picture box
 this.Controls.Add(pb)

 // set memu property
 miFileSave.Index <- 0
 miFileSave.Text <- "&Save"
 miExit.Index <- 1
 miExit.Text <- "&Exit"
 mFile.Index <- 0
 mFile.Text <- "&File"
 mFile.MenuItems.Add(miFileSave) |> ignore
 mFile.MenuItems.Add(miExit) |> ignore
 mainmenu.MenuItems.Add(mFile) |> ignore
 // add menu
 this.Menu <- mainmenu
And add event handler for click event of the save menu item:

// savefile dialog

filesave.DefaultExt <- "pgm"
filesave.Filter <- "PGM Files|*.pgm"
let filesaveClick e =
 if pb.Image = null then
     ()
 else
     if filesave.ShowDialog() = DialogResult.OK then
         let bm = pb.Image :?> Bitmap
         let pic = Matrix.init (bm.Height) (bm.Width) (
             fun i j ->
                 let c = bm.GetPixel(j,i)
                 float c.R * 0.3 + float c.G * 0.59 + float c.B * 0.11
             )
              
         writePgm filesave.FileName pic
miFileSave.Click.Add(filesaveClick)

Put them together:

open System.Windows.Forms
open System.Drawing

type PGMViewer(pic:matrix, name:string) as this =
 inherit Form(Width = 300, Height = 300)

 let mainmenu = new MainMenu()
 let mFile = new MenuItem()
 let miFileSave = new MenuItem()
 let miExit = new MenuItem()
 let pb = new PictureBox()
 let filesave = new SaveFileDialog()

 do
     // set the property of picture box
     let nrow = pic.NumRows
     let ncol = pic.NumCols
     let image = new Bitmap(ncol, nrow)
     printfn "%A %A" nrow ncol
     for i=0 to nrow-1 do
         for j=0 to ncol-1 do
             image.SetPixel(j, i, Color.FromArgb(int pic.[i,j], int pic.[i,j], int pic.[i,j]))
     pb.Image <- image
     pb.Height <- nrow
     pb.Width <- ncol
  
     // add picture box
     this.Controls.Add(pb)
  
     // set memu property
     miFileSave.Index <- 0
     miFileSave.Text <- "&Save"
     miExit.Index <- 1
     miExit.Text <- "&Exit"
     mFile.Index <- 0
     mFile.Text <- "&File"
     mFile.MenuItems.Add(miFileSave) |> ignore
     mFile.MenuItems.Add(miExit) |> ignore
     mainmenu.MenuItems.Add(mFile) |> ignore
     // add menu
     this.Menu <- mainmenu
  
     // savefile dialog
     filesave.DefaultExt <- "pgm"
     filesave.Filter <- "PGM Files|*.pgm"
     let filesaveClick e =
         if pb.Image = null then
             ()
         else
             if filesave.ShowDialog() = DialogResult.OK then
                 let bm = pb.Image :?> Bitmap
                 let pic = Matrix.init (bm.Height) (bm.Width) ( fun i j ->
                         let c = bm.GetPixel(j,i)
                         float c.R * 0.3 + float c.G * 0.59 + float c.B * 0.11
                     )
                      
                 writePgm filesave.FileName pic
     miFileSave.Click.Add(filesaveClick)
  
  
     // set title of the dialog
     this.Text <- "PGMViewer" + " - " + name

And let’s test the viewer now:

(new PGMViewer(readPgm @"C:\temp\temp4.pgm")).Show()

(Copyright: this face is from the ORL face database.)

Attachment

wait.

Matrix and linear algebra in F#, Part II: doing linear algebra via math providers

2010-03-28T23:51:00.001+08:00

Long story short

Suppose we are using the latest release, 1.9.9.9, which does not support linear algebra anymore in its PowerPack. But linear algebra is supported before in 1.9.7.8. Under its bin\ folder, there is a DLL named FSharp.PowerPack.Math.Providers.dll. “Math providers” means F# calls other math libraries, to be specific, an incomplete managed implementation in F#, Netlib’s Lapack+Blas and Intel MKL. The F# one is immature, while the last two are really really good! Plus Netlib’s implementation is free.

However, we cannot just use that provider dll because it has conflicts with latest PowerPack. We need to build a new one! Here are the steps:

1. Compile Netlib’s Lapack to get lapack.dll and blas.dll. Refer to this post of mine. You can also search online to download these two dlls, e.g. http://www.stanford.edu/~vkl/code/libs.html (I haven't tested this).

2. Get the source code of 1.9.7.8. Go to folder FSharp-1.9.7.8\source\fsppack\FSharp.PowerPack. Remove the old reference to FSharp.PowerPack and add two 1.9.9.9 DLLs to the reference:

FSharp.PowerPack.dll and FSharp.PowerPack.Compatibility (Because in 1.9.9.9, Compatibility module is separated out.)

Compile and you get only a few errors regarding “permutation”:

you could just define it explicitly: type permutation = int –> int

Done! Now we move on to

Use Lapack library

First make sure your project

1. references to

FSharp.PowerPack.Compatibility, FSharp.PowerPack and FSharp.PowerPack.Math.Providers.dll.

2. has lapack.dll and blas.dll in the binary output folder. E.g. Debug\ or Release\. (Technically we only need to put them under a folder .Net platform is searchable, similar to Java’s CLASSPATH.)

Finally, let’s write a simple program:

open Microsoft.FSharp.Math
let r = new System.Random()
let A = Matrix.init 4 4 (fun i j -> r.NextDouble() * 1000.0)
let Ainv = Experimental.LinearAlgebra.Inverse A
let eigen = Experimental.LinearAlgebra.EigenValues A

Although this program compiles, it gives a Not-Implemented exception on EigenValues function. (It can compute Inverse of A correctly.) This is because functions in LinearAlgebra module use managed implementation (see linear_algebra_managed.fs in lapack folder) on default. However, the managed implementation is incomplete, and not thoroughly tested.

We need to first open the lapack service:

let isSucc = Experimental.LinearAlgebra.Lapack.Start()

Put this before any computation. (Because there’s no documentation at all, I’ve read all the code in lapack folder to find this line…)

The returned bool value isSucc indicates whether the Lapack service providers successfully loads dlls for Netlib or not.

If you have Intel MKL

Put mkl_def.dll and mkl_lapack.dll under a .Net searchable folder. Notice that the wrapper is for MKL 9.1 not for 10 series.

Btw, the Lapack service provider will choose MKL if both MKL and Netlib are presented.

The Long story

.Net do have some good math libraries. They are efficient (native performance), easy to use. However, they are not free. For the free ones, there’s a Math.Net open source project, which is still under its alpha-testing phase. These open source projects seem to tend to implement all algorithms from starch in C#, rather than link to some mature ones, like Netlib’s lapack. This certainly has some advantages, like more safe code, better memory management in a .Net sense and exception handling. However, this strategy would have a long development time.

I also suspect their efficiency. I have a L-BFGS based logistic regression solver written C++, it runs 3~4 times slower if compiled into .Net using C++/CLI. Drawing a conclusion from this single case is too assertive. (I used STL and STL is slow in C++/CLI, in C++/CLI we’d better use CLI STL.) However, I still trust native libraries to perform basic linear algebra for the performance's seek.

I started to use F# from its 1.9.7.8 (The October CTP, 2009). In that release, the math provider library is already spitted out from the main PowerPack as a standalone DLL. But there’s no documentation on it, except the source code. Other math types, e.g. matrix, vector, are documented. So I didn’t know there exists such a library at first. By occasion, I found the lapack folder in the F# source. Today, I finally have some time to dig into it, get it compiled and know how to start the linear algebra service provider. Fortunately, once started, the service provider is able to automatically find the DLLs of Netlib or MKL, which is really nice.

The code is still compliable using 1.9.9.9, but I am not sure it will in the future as there are a lot of warnings mentioning deprecated language usages, e.g. use OCaml-style ‘^’to connect two strings.

As this code is under MS open source license, there would be license problems if we maintain and distribute it. So a possible solution is to rewrite the service provider and make it open source (e.g. MIT or BSD license), or, persuade MSR Cambridge or MS VS Team to continue the development.

Compiling Lapack for .Net usage

2010-03-28T22:28:00.001+08:00

The aim: we need lapack.dll and blas.dll

A stable and efficient implementation of the common linear algebra routines is essential to engineering. Lapack is the one. Lapack are blas are kind of standard, there are different vendors implementing them, e.g. Intel and AMD both have tuned versions targeted for their hardware.

There are two major implementations available on Windows-Intel platform: Netlib and Intel Math Kernel Library (MKL). Although MKL is very good (Matlab uses it for its basic linear algebra computation), it costs money.Today I’d like to write a post on how to build Netlib’s implementation for .Net usage.

There are two tutorials online for the same purpose:

http://www.codingday.com/compile-lapack-and-blas-as-dll-on-windows

http://icl.cs.utk.edu/lapack-for-windows/lapack/index.html

The first one is quite good. we may encounter some minor problems if we directly follow that guide. I’d like to give a new summary here:

1. Download Lapack+Blas(a single package) from Netlib. Use lapack-3.1.1.tgz.

2. Install MinGw, we only need to install the basic tools + g77 compiler. Because this free compiler only supports Fortran 77 and the latest Lapack is using Fortran 90 now, so in the first step, we use an old version (2007) of Lapack.

Set the PATH for MinGw. And make sure to successfully compile a hello world Fortran program.

3. Copy dlamch.f and slamch.f from INSTALL directory SRC directory

4. Compile BLAS:

g77 --shared -o blas.dll BLAS\SRC\*.f -O3

5. Compile Lapack:

cd SRC
g77 -c *.f -O3
cd ..
g77 -shared -o lapack.dll SRC/*.o blas.dll -O3

If you see some “collect2.exe” error, then change TMP environment variable(originally %USERPROFILE%\AppData\Local\Temp) to a short one, e.g. “c:\temp”

We should have lapack.dll and blas.dll in the Lapack folder now.

Once we have the two native DLLs, we can use P/Invoke in .Net to call linear algebra operations in Lapack. If you use F#, then we don’t even need to know P/Invoke because the F# team already did that for us.

A Note: P/Invoke in Mono platform is also quite convenient, compile lapack and blas into .so files. http://www.mono-project.com/Interop_with_Native_Libraries.

Matrix and linear algebra in F#, Part I: the F# Matrix type

2010-03-28T22:08:00.001+08:00

Every language has libraries, besides the big .Net libraries, F# has two own: the Core, which is shipped with Visual Studio 2010, and the PowerPack, which is an external library developed by MSR Cambridge and Visual Studio Team. Notice that the code quality in PowerPack is actually quite high, it is put outside the Core library because it is evolving fast. Once stable, they may be put into the Core.

Our concern is matrix and linear algebra operations. There is a matrix class in F# PowerPack. However, Microsoft didn’t officially put the documentation online. For F# 1.9.6, there is a outdated page on MSR’s website. But it doe not matter we use the old documentation, since the interface for Math haven’t change much since then.

The namespace for Math is:

Namespace Microsoft.FSharp.Math

In this namespace, we have:

1. complex numbers

2. Big rational numbers

3. vector and row-vector

4. matrix

In this post, I focus on matrix.

The real matrix:

First, Names! There is a type called matrix, which is a matrix holding double or 32-bit long float values.

There is a module called Matrix, inside which there are lots of functions to operate on an F# matrix.

There is also a function/val called matrix, which is used like a constructor to construct a new matrix from lists or arrays.

We can easily create two 3-by-3 matrices using the matrix function:

let A = matrix [ [ 1.0; 7.0; 2.0 ];
               [ 1.0; 3.0; 1.0 ];
               [ 2.0; 9.0; 1.0 ]; ]

let B = matrix [ [ 10.0; 70.0; 20.0 ];
               [ 10.0; 30.0; 10.0 ];
               [ 20.0; 90.0; 10.0 ]; ]

All the member functions and operators associated with matrix type are documented here. Here are some examples:

A+B
A-B
A*B // matrix product
A.*B  // element-wise product
A * 2.0 // scalar product
2.0 * A // this is also ok
-A // negation of a matrix

let b = vector [5.;8.;9.]; // defines a vector
A*b // matrix-vector product

You can get the properties using member functions:

let dim = A.Dimensions
// val dim : int * int = (3, 3), the dimension is a tuple
let A' = A.Transpose // you can use ' in a variable name!
let nrow = A.NumRows
let ncol = A.NumCols
let Anew = A.Copy() // get a new matrix
let Aarr = A.ToArray2D() // convert to a Array2D type
let Avec = A.ToVector() // take the first column of A
let Arvec = A.ToRowVector() // take the fisrt row of A
// ToVector and ToRowVector is usually used
// when you know your matrix is actually a vector

Accessing a matrix

We can have Matlab like access to an F# matrix. One different thing is that the index starts from 0, not 1. Mathematicians like the index to start with 1, e.g. in R and Matlab. While programs like 0-based index, e.g. Numpy for Python.

// notice that the index starts at 0
A.[2,2] //The operator [,] allows access a specific
      //element in the matrix, shorthand for A.Item
A.[2,2] <- 100.0 // change a value
A.[1..2,1..2] // get a sub matrix, shorthand for A.GetSlice
A.[1..2,1..2] <- matrix [[2.;3.]; [8.;9.;]] // set a sub matrix, shorthand for A.SetSlice

We also have 4 member functions: Column, Columns, Row and Rows:

A.Column 2 // Vector<float> = vector [|2.0; 1.0; 1.0|]
A.Row 2 // RowVector<float> = rowvec [|2.0; 9.0; 1.0|]
A.Columns (1,2) // starts at column 1, take 2 columns
//val it : Matrix<float> = matrix [[7.0; 2.0]
//                                 [3.0; 1.0]
//                                 [9.0; 1.0]]
A.Rows (1,2) // starts at row 1, take 2 columns

The Matrix module

Similar to that F# list type has a List module containing handy functions like map, fold and etc, the real matrix type matrix also has a module.

let Asum = Matrix.sum A // sum of all elements in A
let Aprod = Matrix.prod A // product of all elements in A
let C = Matrix.create 10 10 1.0 // create a matrix with 1s
let table = Matrix.init 9 9 (fun i j -> (float i + 1.) * (float j + 1.))
// create a matrix with a function
let I10 = Matrix.identity 10 // 10 1s one diagnal
let Atrace = Matrix.trace A // trace sum
let Asqr = Matrix.map (fun x -> x*x) A // A^2

And there are some repetitions on the member function/operators of matrix type. E.g. Matrix.add, Matrix.set, Matrix.get, Matrix.toVector, Matrix.toRowVector, Matrix.transpose, etc.

Sparse matrix

let D = Matrix.initSparse 3 3 [ (0,0,1.0); (1,1,2.0); (2,2,3.0); ]

// init a sparse 3-by-3 matrix
//val it : matrix = matrix [[1.0; 0.0; 0.0]
//                          [0.0; 2.0; 0.0]
//                          [0.0; 0.0; 3.0]]
D.IsSparse // sparse test

let E = Matrix.initSparse 100000 100000 [ (0,0,1.0); (1,1,2.0); (2,2,3.0); ]

let Esum = Matrix.sum E

However, map, fold, exists, .. are not supported on sparse matrix.

Int Matrix, BigNum Matrix, and others

To know about Generic matrix, you may want to read another post of mine:

http://fdatamining.blogspot.com/2010/03/f-inumerics-interface-and-matrix-class.html

which also discusses some implementation details of the Matrix class, and how to define your own matrix, e.g. a Pixel matrix.

F# INumerics Interface, and Matrix class

2010-03-27T22:31:00.001+08:00

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)

Some Sort Algorithms in F#

2010-03-27T00:51:00.001+08:00

// insertSort

let rec insert x lst =
   match lst with
   | [] -> [x]
   | y::rest -> if x < y then x::y::rest else y::(insert x rest)

let rec insertSort = function
   | [] -> []
   | x::rest -> insert x (insertSort rest)
  
  
// mergeSort
let rec split = function
   | [] -> ([], [])
   | a::[] -> ([a], [])
   | a::b::rest -> let p, q = split rest in (a::p, b::q)

let rec merge p q =
   match p, q with
   | [], [] -> []
   | a, [] -> a
   | [], a -> a
   | x::xs, y::ys when x<=y -> x::merge xs q
   | x::xs, y::ys -> y::merge p ys
  
let rec mergeSort = function
   | [] -> []
   | a::[] -> [a]
   | lst ->
       let p, q = split lst
       let ps = mergeSort p
       let qs = mergeSort q
       merge ps qs
      
// quickSort
let rec quickSort = function
   | [] -> []
   | [a] -> [a]
   | x::xs ->
       quickSort [for y in xs do if y<=x then yield y]
       @ [x]
       @ quickSort [for y in xs do if y>x then yield y]

test from live writer

2010-03-17T14:39:00.001+08:00

seems to be good!

F# online videos

2009-12-13T15:06:00.000+08:00

I found the following videos are very good tutorials for F#.

Luke Hoban
http://microsoftpdc.com/Sessions/FT20
http://blogs.msdn.com/lukeh/

Luca Bolognese
http://channel9.msdn.com/pdc2008/TL11/
http://blogs.msdn.com/lucabol/

Don Syme
http://blogs.msdn.com/dsyme/
Don, the creator of F#, has some videos too:
http://channel9.msdn.com/posts/martinesmann/Don-Syme-FSharp-and-functional-programming-in-NET/

-- some special videos--
Big Algorithms Made Easy with Microsoft's F# by Joel Pobar
http://www.slideshare.net/chlong/dev450-pobar-2884153
http://www.msteched.com/online/view.aspx?tid=3cd60332-dfcb-4b53-b8c5-8aee09648139

InfoQ, Don Syme interview:
http://www.infoq.com/interviews/don-syme-fsharp-2_0

http://www.infoq.com/presentations/Pragmatic-F-Sharp-in-Action