Using Numerical Libraries for .NET from F#
With the release of version 5.0 of our Numerical Libraries for .NET, we’ve also considerably improved our support for the F# language. In this article I want to walk you through setting up F# Interactive to use the Extreme Optimization Numerical Libraries for .NET, and show off some of the new features. This tutorial will focus on F# 3.0 in Visual Studio 2012. Most of the functionality also works with F# 2.0 in Visual Studio 2010.
Setting up F# Interactive
The first step is to bring the Extreme Optimization Numerical Libraries for .NET assemblies into the F# Interactive environment:
#I @”c:\Program Files\Extreme Optimization\Numerical Libraries for .NET\bin\Net40\”;; #r “Extreme.Numerics.Net40.dll”;; #r “Extreme.Numerics.FSharp30.Net40.dll”;;
open Extreme.Mathematics;;
fsi.AddPrinter(fun (A : Extreme.Mathematics.ISummarizable) -> A.Summarize());;
We do this by first adding the path to the Extreme Optimization binaries to the DLL search path using the #I command. (Note: the path depends on the installation folder and may be different on your machine.) We the main assembly Extreme.Numerics.Net40, and the F# specific assembly Extreme.Numerics.FSharp30.Net40. Finally, we tell fsi to use a pretty printer for objects that implement the ISummarizable interface. This will print out vectors and matrices in a user-friendly form.
Working with vectors and matrices
So let’s get our hands dirty right away. When we referenced the Extreme.Numerics.FSharp30.Net40 dll, this automatically opened several modules that make it easier to do linear algebra in F#. We can call Vector and Matrix factory methods, or we can use some convenient F# functions:
> let A = rand 3 3;;
val A : LinearAlgebra.DenseMatrix = 3x3 DenseMatrix [[0.0780,0.1007,0.6917] [0.5422,0.8726,0.5718] [0.3528,0.0631,0.0772]]
let c = !![1.0..10.0];; val c : LinearAlgebra.DenseVector = [1.0000,2.0000,3.0000,4.0000,5.0000,6.0000,7.0000,8.0000,9.0000,10.0000]
Notice the !! operator, which creates a vector from a list or sequence. We can access elements of vectors and matrices using the familiar F# syntax. Slices work too, as do integer sequences and predicates:
> A.[0,1];; val it : float = 0.100694183
A.[1,1] A.[1..2,0..1];; val it : Matrix = 2x2 DenseMatrix [[0.5422,99.0000] [0.3528,0.0631]] c.[1..3];; val it : Vector = [2.0000,3.0000,4.0000] c.[fun x -> sin x > 0.0];; val it : Vector = [1.0000,2.0000,3.0000,7.0000,8.0000,9.0000] c.[[2;4;5]];; val it : Vector = [3.0000,5.0000,6.0000] c.[[for i in 1..7 -> i % 3]];; val it : Vector = [2.0000,3.0000,1.0000,2.0000,3.0000,1.0000,2.0000]
Let’s do some simple calculations:
> A*b;; val it : LinearAlgebra.DenseVector = [3.7378,5.1463,0.8651]
A.T*b;; val it : Vector = [2.9263,2.1614,2.2214] 2.0 * b;; val it : LinearAlgebra.DenseVector = [2.0000,4.0000,10.0000] 1.0 + b;; val it : Vector = [2.0000,3.0000,6.0000] log (1.0 + b);; val it : Vector = [0.6931,1.0986,1.7918] det A;; val it : float = -0.1707208593 inv A;; val it : Matrix = 3x3 DenseMatrix [[-0.1835,-0.2101,3.1982] [-0.9361,1.3939,-1.9355] [1.6027,-0.1792,-0.0791]] let b = Vector.Create(3, fun i -> 1.0 + (float i)**2.0);; val b : LinearAlgebra.DenseVector = [1.0000,2.0000,5.0000] let x = A.Solve(b);; val x : LinearAlgebra.DenseVector = [15.3876,-7.8258,0.8488]
The last command solved the system of equations represented by the matrix A and solved it for the right-hand side b. We can also work with matrix decompositions. There are actually two ways. Let’s use the SVD as an example. First, we can get the factors of the decomposition.
> let (U,Σ,V) = A.svd();;
val Σ : LinearAlgebra.DiagonalMatrix = 3x3 DiagonalMatrix [[1.2980,0.0000,0.0000] [0.0000,0.5101,0.0000] [0.0000,0.0000,0.2578]] val V : Matrix = 3x3 DenseMatrix [[-0.4492,0.4032,-0.7973] [-0.6411,0.4760,0.6020] [-0.6222,-0.7816,-0.0446]] val U : Matrix = 3x3 DenseMatrix [[-0.4083,-0.9041,-0.1260] [-0.8928,0.3667,0.2616] [-0.1903,0.2193,-0.9569]]
Or, if we just want the singular values, we can get those as well:
> let σ = A.svd();;
val σ : LinearAlgebra.DenseVector = [1.2980,0.5101,0.2578]
The second way to work with decompositions is by creating decomposition objects. Let’s work with a non-square matrix this time:
> let C = rand 15 4;;
val C : LinearAlgebra.DenseMatrix = 15x4 DenseMatrix [[0.4198,0.8147,0.8530,0.5364] [0.9181,0.7706,0.3164,0.8773] [0.4768,0.1423,0.7419,0.6521] [0.2463,0.4579,0.3474,0.0311] [0.0432,0.6366,0.4928,0.2399] … [0.3991,0.3426,0.4452,0.1276] [0.7867,0.3247,0.5256,0.1940] [0.6504,0.0943,0.1169,0.4266] [0.0436,0.6716,0.5230,0.4922] [0.5329,0.5422,0.2448,0.0547]]
let d = (rand 15 1).GetColumn(0);;
val d : Vector = [0.7446,0.8430,0.6296,0.1210,0.2190,…,0.8281,0.6643,0.0355,0.7120,0.3582]
Here we see that the output is truncated, and only the 5 first and last rows are shown. We can find the least squares solution to C*y = d. We’ll also compute the residuals and their norm:
> let y = svdC.GetPseudoInverse() * d;;
val y : Vector = [0.3221,0.4458,0.3383,0.1298]
let r = d - C * y;;
val r : Vector = [-0.1120,-0.0171,0.0770,-0.2840,-0.2765,…,0.3796,0.0632,-0.3110,0.1577,-0.1450]
norm r;; val it : float = 0.8214250504
Working with functions
You’d expect a functional language to work with functions. So let’s see what we can do.
The modules discussed in this section are not automatically opened when you load the F# compatibility assembly. The available modules are:
- Special: special functions, such as Gamma functions, Bessel functions, etc.
- Functional: numerical integration and differentiation, finding zeros, maxima and minima of functions.
- Statistics: some functions that are useful when doing statistical analysis.
- Random: random number streams.
These are only the modules with F# specific functions. The full functionality of the libraries is of course available from F# as well.
Here we’ll focus mainly on the second category. Special functions can be very handy, however, so let’s start with those:
> open Extreme.Numerics.FSharp.Special;;
gamma 0.5;; val it : float = 1.772453851 BesselJ 0.0 1.0;; val it : float = 0.7651976866 let J0 = BesselJ 0.0;;
val J0 : (float -> float)
J0 1.0;; val it : float = 0.7651976866
Notice how we used partial application of the Bessel function to get a function that returns the Bessel function of order 0. We can do something similar in many other places, for example with binomial coefficients:
> let binom9 = binomial 9;;
val binom9 : (int -> float)
[for i in 0..9 -> binom9 i];; val it : float list = [1.0; 9.0; 36.0; 84.0; 126.0; 126.0; 84.0; 36.0; 9.0; 1.0]
Now to the functional stuff. Two of the most basic operations on functions are differentiation and integration. Differentiation is done using the d function. Integration is done using integrate1. Partial application is supported, so you can create a function that computes the numerical derivative at any value you give it.
> open Extreme.Numerics.FSharp.Functional;;
let J1 = BesselJ 1.0;;
val J1 : (float -> float)
d J0 0.5;; val it : float = -0.2422684577 J1 0.5;; // d J0 = -J1 val it : float = 0.2422684577 let dJ0 = d J0;;
val dJ0 : (float -> float)
dJ0 0.5;; val it : float = -0.2422684577 dJ0 1.5;; val it : float = -0.5579365079
Integration works in a similar way. The bounds of the integration interval are supplied as a tuple. Infinite intervals are supported as well.
> integrate1 J1 (0.0,5.0);; val it : float = 1.177596771
J0 5.0 - J0 0.0;; // Should give - the previous result: val it : float = -1.177596771 let iJ1 = integrate1 J1;;
val iJ1 : (float * float -> float)
iJ1 (0.0,5.0);; val it : float = 1.177596771 integrate1 (fun x -> exp -x) (0.0,infinity);; val it : float = 1.0
Finding zeros of functions is easy. Let’s quickly find the first 5 zeros of the Bessel function J0:
> [for x0 in [2.0;5.0;8.0;11.0;14.0] -> findzero (fun x -> J0 x) x0];; val it : float list = [2.404825558; 5.52007811; 8.653727913; 11.79153444; 14.93091771]
Finding the minimum or maximum of a function can be difficult because most software requires that you supply the gradient of the function. Although computing derivatives is high school calculus, it’s still very error prone. Automatic differentiation comes to the rescue here. We’ll use the SymbolicMath class, which contains all kinds of functions that take advantage of symbolic computations to obtain a solution.
One of the most famous optimization problems is the Rosenbrock function. We can find the minimum of this function in just one line of code:
> SymbolicMath.FindMinimum((fun x y -> (1.0 - x)**2.0 + 100.0*(y-x**2.0)**2.0), !![-1.0;-1.0]);; val it : SolutionReport = Extreme.Mathematics.SolutionReport`1[Extreme.Mathematics.Vector] {ConvergenceTest = null; Error = 2.778930274e-09; EvaluationsNeeded = 74; IterationsNeeded = 26; Result = [1.0000,1.0000]; Status = Converged;}
It’s time now for you to play. Download the trial version of the Extreme Optimization Numerical Libraries for .NET here and have fun!