Mathias Brandewinder on .NET, F#, VSTO and Excel development, and quantitative analysis / machine learning.
by Mathias 21. February 2015 12:23

A few weeks ago, I came across DiffSharp, an automatic differentiation library in F#. As someone whose calculus skills have always been rather mediocre (thanks Wolfram Alpha!), but who needs to deal with gradients and the like on a regular basis because they are quite useful in machine learning and numerical methods, the project looked pretty interesting: who wouldn’t want exact and efficient calculations of derivatives? So I figured I would take a couple of hours to experiment with the library. This post is by no means an in-depth evaluation, but rather intended as “notes from the road” from someone entirely new to DiffSharp.

Basics

Suppose I want to compute the derivative of f(x) = √ x at, say, 42.0. Double-checking Wolfram Alpha confirms that f has derivative f’(x) = 1 / (2 x √ x) .

Once DiffSharp is installed via Nuget, we can automatically evaluate f’(x) :

#r @"..\packages\DiffSharp.0.5.7\lib\DiffSharp.dll"
open DiffSharp.AD.Forward

let f x = sqrt x
diff f 42. |> printfn "Evaluated: %f"
1. / (2. * sqrt 42.)  |> printfn "Actual: %f"

Evaluated: 0.077152
Actual: 0.077152

val f : x:Dual -> Dual
val it : unit = ()

First off, obviously, it worked. Without any need for us to perform anything, DiffSharp took in our implementation of f, and computed the correct value. This is really nifty.

The piece which is interesting here is the inferred signature of f. If I were to remove the line that immediately follows the function declaration, f would have the following signature:

val f : x:float –> float

The moment you include the line diff f 42., the inferred type changes drastically, and becomes

val f : x:Dual –> Dual

This is pretty interesting. Because we call diff on f, which expects Duals (a type that is defined in DiffSharp), our function isn’t what we originally defined it to be – and calling f 42.0 at that point (for instance) will fail, because 42.0 is a float, and not a Dual. In other words, DiffSharp leverages type inference pretty aggressively, to convert functions into the form it needs to perform its magic.

Edit: Atilim Gunes Baydin suggested another way around that issue, which is inlining f. The following works perfectly well, and allows to both differentiate f, and use this against floats:

let inline f x = sqrt x
let f' = diff f
f 42.

Thanks for the input!

This has a couple of implications. First, if you work in a script, you need to be careful about how you send your code to the F# interactive for execution. If you process the sample code above line by line in FSI, the evaluation will fail, because f will be inferred to be float –> float. Then, you will potentially need to annotate your functions with type hints, to help inference. As an example, the following doesn’t work:

let g x = 3. * x
diff g 42.

As is, g is still inferred to be of type float –> float, because of the presence of the constant term, which is by default inferred as a float. That issue can be addressed at least two ways – by explicitly marking x or 3. as dual in g, like this:

let g x = (dual 3.) * x
let h (x:Dual) = 3. * x

That’s how far we will go on this – if you want to dive in deeper, the Type Inference page discusses the topic in much greater detail

A tiny example

So why is this interesting? As I mentioned earlier, differentiation is used heavily in numeric algorithms to identify values that minimize a function, a prime example being the gradient descent algorithm. The simplest example possible would be finding a (local) minimum of a single-argument function: starting from an arbitrary value x, we can iteratively follow the direction of steepest descent, until no significant change is observed.

Here is a quick-and-dirty implementation, using DiffSharp:

let minimize f x0 alpha epsilon =
    let rec search x =
        let fx' = diff f x
        if fx' < epsilon
        then x
        else
            let x = x - alpha * fx'
            search x
    search x0

Because DiffSharp handles the differentiation part automatically for us, with only 10 lines of code, we can now pass in arbitrary functions we want to minimize, and (with a few caveats…), and get a local minimum, no calculus needed:

let epsilon = 0.000001

let g (x:Dual) = 3. * pown x 2 + 2. * x + 1.
minimize g 0. 0.1 epsilon |> printfn "Min of g at x = %f"

let h (x:Dual) = sin (pown (x + 30.) 2) 
minimize h 0. 0.1 epsilon |> printfn "Min of h at x = %f"

Min of g at x = -0.333333
Min of h at x = -0.397480

Let’s make sure this is reasonable. g is a quadratic function, which has a minimum or maximum at –b/2*a, that is, –2 / 2 x 3 - this checks out. As for h, inspecting the function plot confirms that it has a minimum around the identified value:

function-plot

Conclusion

I have barely started to scratch the surface of DiffSharp in this post, but so far, I really, really like its promise. While I limited my examples to single-variable functions, DiffSharp supports multivariate functions, and vector operations as well. The way it uses type inference is a bit challenging at first, but seems a reasonable price to pay for the resulting magic. My next step will probably be a less tiny example, perhaps a logistic regression against realistic data. I am very curious to try out the algebra bits – and also wondering in the back of my head how to best use the library in general. For instance, how easy is it to construct a function from external data, and turn it into the appropriate types for DiffSharp to work its magic? How well does this integrate with other libraries, say, Math.NET? We’ll see!

In the meanwhile, I’d recommend checking out the project page, which happens to also be beautifully documented! And, as always, you can ping me on twitter for comments or question.

by Mathias 11. January 2015 18:37

I had the great pleasure to speak at CodeMash this week, and, on my way back, ended up spending a couple of hours at the Atlanta airport waiting for my connecting flight back to the warmer climate of San Francisco – a perfect opportunity for some light-hearted coding fun. A couple of days earlier, I came across this really nice tweet, rendering the results of an L-system:

I ended up looking up L-systems on Wikipedia, and thought this would make for some fun coding exercise. In a nutshell, a L-system is a grammar. It starts with an alphabet of symbols, and a set of rules which govern how each symbol can be transformed into another chain of symbols. By applying these rules to a starting state (the initial axiom), one can evolve it into a succession of states, which can be seen as the growth of an organism. And by mapping each symbol to operations in a logo/turtle like language, each generation can then be rendered as a graphic.

So how could we go about coding this in F#? If you are impatient, you can find the final result as a gist here.

First, I started with representing the core elements of an L-System with a couple of types:

type Symbol = | Sym of char

type State = Symbol list

type Rules = Map<Symbol,State>

type LSystem = 
    { Axiom:State
      Rules:Rules }

A symbol is a char, wrapped in a single-case discriminated union, and a State is simply a list of Symbols. We define the Rules that govern the transformation of Symbols by a Map, which associates a particular Symbol with a State, and an L-System is then an Axiom (the initial State), with a collection of Rules.

Let’s illustrate this on the second example from the Wikipedia page, the Pythagoras tree. Our grammar contains 4 symbols, 0, 1, [ and ], we start with a 0, and we have 2 rules, (1 → 11), and (0 → 1[0]0). This can be encoded in a straightforward manner in our domain, like this:

let lSystem =
    { Axiom = [ Sym('0') ]
      Rules = [ Sym('1'), [ Sym('1'); Sym('1') ]
                Sym('0'), [ Sym('1'); Sym('['); Sym('0'); Sym(']'); Sym('0') ]]
              |> Map.ofList }

Growing the organism by applying the rules is fairly straightforward: given a State, we traverse the list of Symbols, look up for each of them if there is a matching rule, and perform a substitution if it is found, leaving it unchanged otherwise:

(*
Growing from the original axiom
by applying the rules
*)

let applyRules (rs:Rules) (s:Symbol) =
    match (rs.TryFind s) with
    | None -> [s]
    | Some(x) -> x

let evolve (rs:Rules) (s:State) =
    [ for sym in s do yield! (applyRules rs sym) ]

let forward (g:LSystem) =
    let init = g.Axiom
    let gen = evolve g.Rules
    init |> Seq.unfold (fun state -> Some(state, gen state))

// compute nth generation of lSystem
let generation gen lSystem =
    lSystem
    |> forward 
    |> Seq.nth gen
    |> Seq.toList

What does this give us on the Pythagoras Tree?

> lSystem |> generation 1;;
val it : Symbol list = [Sym '1'; Sym '['; Sym '0'; Sym ']'; Sym '0']

Nice and crisp – that part is done. Next up, rendering. The idea here is that for each Symbol in a State, we will perform a substitution with a sequence of instructions, either a Move, drawing a line of a certain length, or a Turn of a certain Angle. We will also have a Stack, where we can Push or Pop the current position of the Turtle, so that we can for instance store the current position and direction on the stack, perform a couple of moves with a Push, and then return to the previous position by a Pop, which will reset the turtle to the previous position. Again, that lends itself to a very natural model:

(*
Modelling the Turtle/Logo instructions
*)

type Length = | Len of float
type Angle = | Deg of float

// override operator later
let add (a1:Angle) (a2:Angle) =
    let d1 = match a1 with Deg(x) -> x
    let d2 = match a2 with Deg(x) -> x
    Deg(d1+d2)

type Inst =
    | Move of Length
    | Turn of Angle
    | Push
    | Pop

let Fwd x = Move(Len(x))
let Lft x = Turn(Deg(x))
let Rgt x = Turn(Deg(-x))

We can now transform our L-system state into a list of instructions, and convert them into a sequence of Operations, in that case Drawing lines between 2 points:

type Pos = { X:float; Y:float; }
type Dir = { L:Length; A:Angle }

type Turtle = { Pos:Pos; Dir:Dir }
type ProgState = { Curr:Turtle; Stack:Turtle list }

let turn angle turtle = 
    let a = turtle.Dir.A |> add angle
    { turtle with Dir = { turtle.Dir with A = a } }

type Translation = Map<Symbol,Inst list>

type Ops = | Draw of Pos * Pos

let pi = System.Math.PI

let line (pos:Pos) (len:Length) (ang:Angle) =
    let l = match len with | Len(l) -> l
    let a = match ang with | Deg(a) -> (a * pi / 180.)
    { X = pos.X + l * cos a ; Y = pos.Y + l * sin a }

let execute (inst:Inst) (state:ProgState) =
    match inst with
    | Push -> None, { state with Stack = state.Curr :: state.Stack }
    | Pop -> 
        let head::tail = state.Stack // assumes more Push than Pop
        None, { state with Curr = head; Stack = tail }
    | Turn(angle) -> 
        None, { state with Curr =  state.Curr |> turn angle }
    | Move(len) -> 
        let startPoint = state.Curr.Pos
        let endPoint = line startPoint len state.Curr.Dir.A
        Some(Draw(startPoint,endPoint)), { state with Curr = { state.Curr with Pos = endPoint } }

let toTurtle (T:Translation) (xs:Symbol list) =

    let startPos = { X = 400.; Y = 400. }
    let startDir = { L = Len(0.); A = Deg(0.) }
    let init = 
        { Curr = { Pos = startPos; Dir = startDir }
          Stack = [] }
    xs 
    |> List.map (fun sym -> T.[sym]) 
    |> List.concat
    |> Seq.scan (fun (op,state) inst -> execute inst state) (None,init)
    |> Seq.map fst
    |> Seq.choose id

We simply map each Symbol to a List of instructions, transform the list of symbols into a list of instructions, and maintain at each step the current position and direction, as well as a Stack (represented as a list) of positions and directions. How does it play out on our Pythagoras Tree? First, we define the mapping from Symbols to Instructions:

let l = 1.
let T = 
    [ Sym('0'), [ Fwd l; ]
      Sym('1'), [ Fwd l; ]
      Sym('['), [ Push; Lft 45.; ]
      Sym(']'), [ Pop; Rgt 45.; ] ]
    |> Map.ofList

… and we simply send that toTurtle, which produces a list of Draw instructions:

> lSystem |> generation 1 |> toTurtle T;;
val it : seq<Ops> =
  seq
    [Draw ({X = 400.0;
            Y = 400.0;},{X = 401.0;
                         Y = 400.0;}); Draw ({X = 401.0;
                                              Y = 400.0;},{X = 401.7071068;
                                                           Y = 400.7071068;});
     Draw ({X = 401.0;
            Y = 400.0;},{X = 401.7071068;
                         Y = 399.2928932;})]

Last step – some pretty pictures. We’ll simply generate a html document, rendering the image using SVG, by creating one SVG line per Draw instruction:

let header = """
<!DOCTYPE html>
<html>
<body>
<svg height="800" width="800">"""

let footer = """
</svg>
</body>
</html>
"""

let toSvg (ops:Ops seq) =
    let asString (op:Ops) = 
        match op with
        | Draw(p1,p2) -> sprintf """<line x1="%f" y1="%f" x2="%f" y2="%f" style="stroke:rgb(0,0,0);stroke-width:1" />""" p1.X p1.Y p2.X p2.Y 

    [ yield header
      for op in ops -> asString op
      yield footer ]
    |> String.concat "\n"

open System.IO

let path = "C:/users/mathias/desktop/lsystem.html"
let save template = File.WriteAllText(path,template)

And we are pretty much done:

> lSystem |> generation 8 |> toTurtle T |> toSvg |> save;;
val it : unit = ()

… which produces the following graphic:

image

Pretty neat! Just for fun, I replicated the Sierpinski Triangle example as well:

let sierpinski () =

    let lSystem =
        { Axiom = [ Sym('A') ]
          Rules = [ Sym('A'), [ Sym('B'); Sym('>'); Sym('A'); Sym('>'); Sym('B') ]
                    Sym('B'), [ Sym('A'); Sym('<'); Sym('B'); Sym('<'); Sym('A') ]]
                  |> Map.ofList }

    let l = 1.
    let T = 
        [ Sym('A'), [ Fwd l; ]
          Sym('B'), [ Fwd l; ]
          Sym('>'), [ Lft 60.; ]
          Sym('<'), [ Rgt 60.; ] ]
        |> Map.ofList

    lSystem 
    |> generation 9
    |> toTurtle T
    |> toSvg 
    |> save

… which results in the following picture:

image

That’s it for tonight! I had a lot of fun coding this (it certainly made the flight less boring), and found the idea of converting code to turtle instructions, with a stack, pretty interesting. Hope you enjoyed it, and if you end up playing with this, share your creations on Twitter and ping me at @brandewinder!

Gist for the whole code here

by Mathias 31. December 2014 10:59

Well, we are in the last hours of 2014, and I am nearly recovered from the craziness that was the F# Europa Tour 2014, so here we go – the Tour, in cold, hard facts (after all, I am a numbers’ guy):

  • 40 days of travelling across Europe.
  • 16 talks.
  • 5 workshops (about 50 hours total).
  • 9 countries.
  • 6991 miles (11,250 kilometers) travelled, roughly (this is straight-line city to city, so the actual number is probably a good deal larger).
  • 14 hours of bus.
  • roughly 50 hours of train.
  • roughly 28 hours of plane.
  • 12 cities visited (and spoken at!).
  • I lost track of how many gallons of beer were ingested. This is big data.
  • 500 attendees? Maybe more? See previous data point.
  • Delivered hundreds of shiny fsharp.org stickers to F# Communities across Europe. [btw, in case you didn't hear - the F# Software Foundation is now a full-fledged, legally established entity, and YOU can be a member. Check it out!]

Now for the important qualitative questions:

  • Where did I eat the best bacon? This came as a surprise to me, but I have to say, the bacon I ate in Dublin, Ireland was amazing. Twice.
  • Where does one find the best beer in Europe? This is a hard one – I had a chance to sample great beers from all over the place. I would say, Munich and its Biergarten rules, but the live beers at BuildStuff in Vilnius, Lithuania, were a very nice surprise.
  • What’s the weirdest thing I ate? This one goes to Norway and its Lutefisk, a traditional Christmas fish dish. It’s definitely a regional specialty, as in, a specialty which didn’t expand beyond a limited regional area, for good reasons. For the record, I actually enjoyed it!
  • What was the worst travelling mistake? Booking a last minute train ticket from Paris to Aarhus, Denmark, to realize in the train that instead of a nice sleeping car, I would be spending 22 hours sitting in a train with no food on board.
  • Biggest scare: every person who has given a talk will tell you, relying on the internet and anything live in a presentation is a rookie mistake. This is great advice, which is why I completely ignored it. It all worked just fine, but learning that Azure had been down for a couple of hours, right before a talk at BuildStuff which 100% required a live deployment to Azure to work, did give me some cold sweat.

Would I do it again? In a heartbeat! It was a bit crazy, and definitely exhausting, but a ton of fun. All of you who helped out making this happen, from the bottom of my heart, thank you! The F# Community is absolutely fantastic, packed with energy and a good, friendly vibe, and everywhere I went felt like family. You all kept me going, so again, thank you (you know who you are)! In the meanwhile, I wish you all a happy year 2015 ahead, let’s make that one even better than 2014, and I hope to see many of you again this year! And, as always, feel free to ping me on Twitter as @brandewinder.

by Mathias 1. March 2014 14:32

During some recent meanderings through the confines of the internet, I ended up discovering the Winnow Algorithm. The simplicity of the approach intrigued me, so I thought it would be interesting to try and implement it in F# and see how well it worked.

The purpose of the algorithm is to train a binary classifier, based on binary features. In other words, the goal is to predict one of two states, using a collection of features which are all binary. The prediction model assigns weights to each feature; to predict the state of an observation, it checks all the features that are “active” (true), and sums up the weights assigned to these features. If the total is above a certain threshold, the result is true, otherwise it’s false. Dead simple – and so is the corresponding F# code:

type Observation = bool []
type Label = bool
type Example = Label * Observation
type Weights = float []

let predict (theta:float) (w:Weights) (obs:Observation) = 
    (obs,w) ||> Seq.zip 
    |> Seq.filter fst 
    |> Seq.sumBy snd 
    |> ((<) theta)

We create some type aliases for convenience, and write a predict function which takes in theta (the threshold), weights and and observation; we zip together the features and the weights, exclude the pairs where the feature is not active, sum the weights, check whether the threshold is lower that the total, and we are done.

In a nutshell, the learning process feeds examples (observations with known label), and progressively updates the weights when the model makes mistakes. If the current model predicts the output correctly, don’t change anything. If it predicts true but should predict false, it is over-shooting, so weights that were used in the prediction (i.e. the weights attached to active features) are reduced. Conversely, if the prediction is false but the correct result should be true, the active features are not used enough to reach the threshold, so they should be bumped up.

And that’s pretty much it – the algorithm starts with arbitrary initial weights of 1 for every feature, and either doubles or halves them based on the mistakes. Again, the F# implementation is completely straightforward. The weights update can be written as follows:

let update (theta:float) (alpha:float) (w:Weights) (ex:Example) =
    let real,obs = ex
    match (real,predict theta w obs) with
    | (true,false) -> w |> Array.mapi (fun i x -> if obs.[i] then alpha * x else x)
    | (false,true) -> w |> Array.mapi (fun i x -> if obs.[i] then x / alpha else x)
    | _ -> w

Let’s check that the update mechanism works:

> update 0.5 2. [|1.;1.;|] (false,[|false;true;|]);;
val it : float [] = [|1.0; 0.5|]

The threshold is 0.5, the adjustment multiplier is 2, and each feature is currently weighted at 1. The state of our example is [| false; true; |], so only the second feature is active, which means that the predicted value will be 1. (the weight of that feature). This is above the threshold 0.5, so the predicted value is true. However, because the correct value attached to that example is false, our prediction is incorrect, and the weight of the second feature is reduced, while the first one, which was not active, remains unchanged.

Let’s wrap this up in a convenience function which will learn from a sequence of examples, and give us directly a function that will classify observations:

let learn (theta:float) (alpha:float) (fs:int) (xs:Example seq) =
    let updater = update theta alpha
    let w0 = [| for f in 1 .. fs -> 1. |]    
    let w = Seq.fold (fun w x -> updater w x) w0 xs
    fun (obs:Observation) -> predict theta w obs

We pass in the number of features, fs, to initialize the weights at the correct size, and use a fold to update the weights for each example in the sequence. Finally, we create and return a function that, given an observation, will predict the label, based on the weights we just learnt.

And that’s it – in 20 lines of code, we are done, the Winnow is implemented.

More...

by Mathias 15. February 2014 12:51

My favorite column in MSDN Magazine is Test Run; it was originally focused on testing, but the author, James McCaffrey, has been focusing lately on topics revolving around numeric optimization and machine learning, presenting a variety of methods and approaches. I quite enjoy his work, with one minor gripe –his examples are all coded in C#, which in my opinion is really too bad, because the algorithms would gain much clarity if written in F# instead.

Back in June 2013, he published a piece on Amoeba Method Optimization using C#. I hadn’t seen that approach before, and found it intriguing. I also found the C# code a bit too hairy for my feeble brain to follow, so I decided to rewrite it in F#.

In a nutshell, the Amoeba approach is a heuristic to find the minimum of a function. Its proper respectable name is the Nelder-Nead method. The reason it is also called the Amoeba method is because of the way the algorithm works: in its simple form, it starts from a triangle, the “Amoeba”; at each step, the Amoeba “probes” the value of 3 points in its neighborhood, and moves based on how much better the new points are. As a result, the triangle is iteratively updated, and behaves a bit like an Amoeba moving on a surface.

Before going into the actual details of the algorithm, here is how my final result looks like. You can find the entire code here on GitHub, with some usage examples in the Sample.fsx script file. Let’s demo the code in action: in a script file, we load the Amoeba code, and use the same function the article does, the Rosenbrock function. We transform the function a bit, so that it takes a Point (an alias for an Array of floats, essentially a vector) as an input, and pass it to the solve function, with the domain where we want to search, in that case, [ –10.0; 10.0 ] for both x and y:

#load "Amoeba.fs"

open Amoeba
open Amoeba.Solver

let g (x:float) y =
    100. * pown (y - x * x) 2 + pown (1. - x) 2

let testFunction (x:Point) =
    g x.[0] x.[1]

solve Default [| (-10.,10.); (-10.,10.) |] testFunction 1000

Running this in the F# interactive window should produce the following:

val it : Solution = (0.0, [|1.0; 1.0|])
>

The algorithm properly identified that the minimum is 0, for a value of x = 1.0 and y = 1.0. Note that results may vary: this is a heuristic, which starts with a random initial amoeba, so each run could produce slightly different results, and might at times epically fail.

More...

Comments

Comment RSS