Tutorial: Using Math.NET Numerics in F#

Applies to: Functional Programming

Authors: Yin Zhu

Get this book in Print, PDF, ePub and Kindle at manning.com. Use code “MSDN37b” to save 37%.

Summary: Math.NET Numerics aims to be the standard open source math library for .NET Framework. It is released under the MIT license and so it can be used in both commercial and noncommercial projects, as well as modified and redistributed. This tutorial gives an introduction on how to use the various components of this library.

This topic contains the following sections.

• Introducing Math.NET Numerics
• Using Vectors and Matrices
• Using Linear Algebra Routines
• Using Probability and Statistics APIs
• Summary
• Additional Resources
• See Also

This article is associated with Real World Functional Programming: With Examples in F# and C# by Tomas Petricek with Jon Skeet from Manning Publications (ISBN 9781933988924, copyright Manning Publications 2009, all rights reserved). No part of these chapters may be reproduced, stored in a retrieval system, or transmitted in any form or by any means—electronic, electrostatic, mechanical, photocopying, recording, or otherwise—without the prior written permission of the publisher, except in the case of brief quotations embodied in critical articles or reviews.

Introducing Math.NET Numerics

Math.NET Numerics is part of a larger project called Math.NET, which is an umbrella project for open source math libraries of the .NET Framework. In the future, the Math.NET project may contain various aspects of computer mathematics, such as symbolic computing, signal processing, and machine learning.

Currently, Math.NET Numerics is the most actively developed part of Math.NET. Among others, it contains a wide range of numerical algorithms including:

• Types representing dense as well sparse vectors and matrices.

• A library containing a full set of standard linear algebra routines.

• A range of random number generators with different strengths.

• Additional features that include support for certain numerical integration and statistical functions.

The Math.NET Numerics library is implemented in the safe subset of C# language, which makes it very portable. The library runs on a wide range of platforms, including Microsoft .NET Framework (for developing web, server, and desktop applications), Silverlight (for developing rich client-side and Windows Phone 7 applications), and Mono as well as Moonlight (for developing cross-platform applications).

Another important feature of Math.NET Numerics is that it provides wrappers for highly-optimized native math libraries such as Intel MKL and AMD ACML. This makes it possible to utilize specialized CPU instruction sets when they are available. For example, the matrix multiplication in Intel MKL runs several times faster than portable C# code. The rest of this article introduces several Math.NET Numerics features in detail.

Using Vectors and Matrices

The Math.NET Numerics library can be downloaded from the official homepage at CodePlex. The two assemblies that need to be referenced when using the library from F# are MathNet.Numerics.dll and MathNet.Numerics.FSharp.dll. When creating a compiled project, the assemblies need to be referenced using the Add Reference dialog in Visual Studio. The download contains the compiled versions of the assemblies for .NET Framework 4.0 and Sliverlight. When using Math.NET Numerics on other platforms, the library may need to be compiled from the source.

The F# support for interactive scripting makes it more convenient to write and explore mathematical problems using an F# Script File and F# Interactive in Visual Studio. Math.NET Numerics can be referenced from an F# Script File using the following commands:

#r @"C:\...\MathNet.Numerics\lib\Net40\MathNet.Numerics.FSharp.dll"
#r @"C:\...\MathNet.Numerics\lib\Net40\MathNet.Numerics.dll"

The MathNet.Numerics.dll assembly contains types for vectors and matrices, linear algebra routines, and others. The MathNet.Numerics.FSharp.dll assembly is a small F# wrapper for Math.NET Numerics, which helps to create vectors and matrices using basic F# types such as sequences and lists. The distribution also comes with documentation (MathNet.Numerics.chm) and a comprehensive set of examples in C#.

Creating and Accessing Vectors and Matrices

The most fundamental data structures of a numerical library are types representing vectors and matrices. In Math.NET Numerics, there are four kinds of vectors and matrices for different numeric types (float, double, complex32, and 64-bit complex). The types are located in four sub-namespaces of the MathNet.Numerics.LinearAlgebra namespace. A script file typically imports vectors and matrices from one of the four sub-namespaces (Single, Double, Complex32 or Complex). The following snippet shows how to use vectors from the Double namespace. For historical reasons, F# calls the double type float, and the single type is called float32:

open MathNet.Numerics.LinearAlgebra.Double
let vector1 = new DenseVector(5)
let vector2 = new DenseVector(5, 3.0)
let vector3 = new DenseVector([| 1.0; 2.0; 3.0; 4.0; 5.0 |])
let vector4 = new DenseVector(vector3)

The DenseVector type can be created using an overloaded constructor. When provided with an integer, it creates a zero vector of the specified size. The second example creates a vector of the given size with each element set to 3.0. Finally, the last two examples create a vector from the data contained in an array or another vector, respectively.

A vector is a mutable type, so it is possible to change its values. The following snippet shows several examples of reading and writing vector elements:

> vector1.[0] <- 100.0;;
> vector1.[0];;
val it : float = 100.0

> let array1 = vector1.ToArray();;
> array1.[0] <- 10.0;;
> vector1.[0];;
val it : float = 100.0

> [ for i, val in vector1.GetIndexedEnumerator() -> (i, val) ];;
val it : (int * float) list =
  [(0, 100.0); (1, 0.0); (2, 0.0); (3, 0.0); (4, 0.0)]
> let subvector = vector1.SubVector(1, 3);;
> subvector.[0] <- 42.0;;
> vector1.[1];; 
val it : float = 0.0

The first two commands demonstrate how to access individual elements using an indexer. As already noted, vector is mutable and the value can be changed using the <- operator. The next couple of commands demonstrate how to turn a vector into a standard F# type. The ToArray method copies the data from the vector to an array. The GetIndexedEnumerator member can be used to get elements together with their indices. Finally, the last few commands demonstrate how to get a subvector from a vector. A subvector is a copy of the original vector, that is, changing the values of the subvector does not affect the original vector.

Similarly to the DenseVector type, there are also four different DenseMatrix types defined in the same four namespaces. The following listing shows how to create double matrices using the constructors of DenseMatrix.

open MathNet.Numerics.LinearAlgebra.Double
let matrix1 = new DenseMatrix(3)
let matrix2 = new DenseMatrix(2, 3)
let matrix3 = new DenseMatrix(2, 3, 3.0)

let arr2d = Array2D.init 2 3 (fun i j -> 
    float (i * 10 + j))
let matrix4 = new DenseMatrix(arr2d)

let arr1d = [| 1.0; 4.0; 2.0; 5.0; 3.0; 6.0 |]
let matrix5 = new DenseMatrix(2, 3, arr1d)

let matrixI = DenseMatrix.Identity(5)

The first three examples create a matrix filled with a constant value. When given a single integer as an argument, DenseMatrix creates a square matrix, so matrix1 is a 3 x 3 matrix containing zeros. The matrix2 value is a 2x3 matrix containing zeros and matrix3 is a 2x3 matrix filled with the value 3.0.

The next two examples create a matrix from an F# array. The matrix4 value is created from a two-dimensional array created using the Array2D.init function. The DenseMatrix constructor can also load data from a single-dimensional array. It assumes that the data is stored in column-major order. Finally, the last example uses the static Identity method. The example creates an identity matrix (containing ones on the diagonal and zeros everywhere else) of size 5 × 5.

Similarly to vectors, matrices are mutable types and the elements can be accessed or modified using the indexer and methods of the Matrix type:

let matrix = new DenseMatrix(3, 3)
for i = 0 to matrix.RowCount - 1 do
    for j = 0 to matrix.ColumnCount - 1 do
        matrix.[i, j] <- float (i * 10 + j)

let sub = matrix.SubMatrix(0, 1, 1, 1)
let row = matrix.Row(2, 1, 2);

The first part of the snippet creates a new matrix filled with zeros and then initializes it using an indexer. The second part of the snippet creates a submatrix of size 2 × 1. The return type of SubMatrix is always Matrix<'T> even if the slice actually represents a row. Matrix<'T> is a base class of DenseMatrix, and it provides a common interface for both sparse and dense matrices. A vector representing a single row or column of the matrix can be obtained using Row and Column members. Note that methods like SubMatrix and Row return new copies of the subpart of the original matrix.

As this section demonstrated, working with vectors and matrices using the standard .NET programming API is quite nice, but the library also provides several extensions for F#.

Using F# Extensions

Math.NET Numerics comes with a small assembly MathNet.Numerics.FSharp.dll that contains several extensions designed especially for F#. Functions from this assembly make creating vectors and matrices in F# even easier. After referencing the assembly, F# functionality can be found in the MathNet.Numerics.FSharp namespace:

> #r @"C:\...\MathNet.Numerics\lib\Net40\MathNet.Numerics.FSharp.dll"
  open MathNet.Numerics.FSharp
  open MathNet.Numerics.LinearAlgebra
  open MathNet.Numerics.LinearAlgebra.Double;;
(...)

> let vec1 = vector [ for i in 0 .. 99 -> float i / 100.0];;
val vec1 : Generic.Vector<float> = seq [0.0; 0.01; 0.02; 0.03; ...]
> let vec2 = DenseVector.init 100 (fun i -> float i / 100.0);;
val vec2 : DenseVector = seq [0.0; 0.01; 0.02; 0.03; ...] 

> let vec3 = 
      seq { for i in 1 .. 100 do yield float i }
      |> DenseVector.ofSeq;;
val vec3 : DenseVector = seq [1.0; 2.0; 3.0; 4.0; ...]

The extension provides a function named vector that can be used to construct vectors using the F# sequence comprehension syntax. Similarly to the standard Array and Array2D modules, the extension also provides a function DenseVector.init that takes a desired vector length and initializes its values using a provided function. Finally, it is also possible to use DenseVector.ofSeq to turn any collection that implements the IEnumerable<'T> type (called seq<'T> in F#) into a vector.

The listings so far demonstrated how to create matrices and vectors and how to access their values. The following section looks at standard operations of linear algebra.

Using Linear Algebra Routines

Linear algebra is at the core of applied mathematics. Linear algebra applications include machine learning, data mining, and artificial intelligence, which are then used to build high-level applications such as character recognition systems and road routing systems. The following sections show how to use linear algebra in Math.NET Numerics and the subsequent sections provide more information about the internals and about using more efficient native implementations.

Writing Linear Algebra Calculations

Working with linear algebra in the Math.NET Numerics library is very easy. Most of the important operations are available as members of the matrixtype. The following code example shows how to create a random matrix and calculate its inverse:

let rnd = new System.Random()
let a = DenseMatrix.init 500 500 (fun _ _ -> rnd.NextDouble())
let a' = a.Inverse()

An inverse matrix can be calculated using the Inverse member. Other linear algebra operations can be accessed using the same syntax and they can be easily explored using Visual Studio's IntelliSense. The next three snippets show how to decompose a matrix using the QR decomposition, Eigenvalue Decomposition (EVD), and Singular Value Decomposition (SVD) algorithms:

let qr = A.QR()
let q, r = qr.Q, qr.R
 
let eigen = A.Evd()
let evals, evecs = eigen.EigenValues(), eigen.EigenVectors()
 
let svd = A.Svd(true)
let u, s, t = svd.U(), svd.S(), svd.VT()

The first two lines show how to calculate the QR decomposition on a dense matrix A. QR decomposition is often used for solving linear least squares problems. The following two lines show how to perform the Eigen decomposition of a matrix, which has many applications, such as face recognition and Principal Components Analysis (PCA). The final two lines show how to calculate the SVD.

Understanding Math.NET Numerics Providers

Math.NET Numerics uses an extensibility model that can be used to provide an alternative implementation of basic linear algebra operations. The extensibility model is based on the concept of the linear algebra provider. Technically, a provider is a .NET interface that is similar to the standard linear algebra operations defined in LAPAK. The following F# code snippet shows a part of the interface:

type ILinearAlgebraProvider<'T> = 
    abstract AddVectorToScaledVector : 
        ('T[] * 'T * 'T[] * 'T[] ) -> unit
    abstract DotProduct : 
        ('T[] * 'T[]) -> unit
    abstract MatrixMultiply: 
        ('T[] * int * int * 'T[] * int * int * 'T[]) -> unit
    abstract SingularValueDecomposition:
        (bool * 'T[] * int * int * 'T[] * 'T[] * 'T[] ) -> unit
    // (Many other functions omitted)

The above listing shows that the style of parameters is very similar to that of LAPAK library. For example, SingularValueDecomposition takes seven parameters, the first bool variable indicates computing the singular U and VT vectors or not, the following 1-D array and two integers represent the input matrix with two integers as the dimensions of the matrix, and the rest three arrays store the decomposition result. However, Math.NET Numerics has a full C# implementation of the above interface. As with the vector and matrix type, there is a C# implementation of this interface for four types: float, double, complex32 and complex.

The design based on providers allows multiple implementations. The Math.NET Numerics library contains several standard implementations that can be found in the MathNet.Numerics.Algorithms.LinearAlgebra namespace. As the previous section demonstrated, the provider is not used directly. All functionality is accessed using members of the Matrix<'T> type.

The default managed implementation is reasonably efficient and works in any environment. When performance is critical, it is possible to use a native implementation such as Intel MKL or ACM ACML. The following section provides more information.

Using Native Algebra Libraries

Most of the time, numerical computations are executed on a desktop or server computer with full rights, so it is possible to call native code. In that case, using a highly optimized native provider for Intel MKL or AMD ACML will give better performance.

Because of license restrictions, Math.NET Numerics cannot provide a precompiled Intel MKL and AMD ACML runtime. This section discusses how to compile and use the AMD ACML library to perform linear algebra. The AMD ACML library is freely downloadable and is free for non-commercial use on Windows. Intel MKL has more restrictions and is available only under a commercial license. The following steps describe how to obtain and compile the AMD ACML library.

Compiling and Using AMD ACML in Math.NET Numerics

1. The source code of the ACML library can be found on the official AMD downloads website. At the time of this writing, the latest version is ACML 4.4. ACML is available on several operating systems compiled with different FORTRAN compilers. This tutorial uses a 32-bit Windows version compiled with Intel FORTRAN compiler. Following the default setting, the library installs to C:\AMD\acml4.4.0.

2. To compile the wrapper, locate the Math.NET Numerics installation package, open the Visual Studio 2010 solution NativeWrappers.sln in src\NativeWrappers\Windows, and compile ACMWrapper project. If the ACML library isn't installed in the default folder, the C++ project properties need to be changed to specify paths to the ACML headers and static libraries.

3. After the complication succeeds, the resulting binary named MathNET.Numerics.ACML.dll can be found in the Release or Debug subdirectories. Copy this assembly to the folder containing MathNET.Numerics.dll so that it can be loaded by Math.NET Numerics.

4. Next, copy libacml_mp_dll.dll, libifcoremd.dll, libiomp5md.dll and libmmd.dll from the ACML installation to same folder. This way, the managed wrapper will be able to locate the native implementation.

After following the above steps, the Math.NET Numerics library can be configured to use the newly compiled provider. This can be done using the following snippet:

open MathNet.Numerics
open MathNet.Numerics.Algorithms.LinearAlgebra

let acml = new Acml.AcmlLinearAlgebraProvider()
Control.LinearAlgebraProvider <- acml

The performance of the two implementations can be compared by measuring the time needed to execute the examples from the section Writing Linear Algebra Calculations. The easiest way to do that is to use the #time directive in F# Interactive:

> #time;;
(...)

> // Example source code goes here 
  ;;

The performance results may vary depending on the configuration. The following numbers are based on measurements made on the author's ThinkPad T61 Laptop (with 2.0 GHz Core Dual CPU and 3.0 GB of memory). Both of the implementations were executed using just a single thread. The average running time of the function (over 10 repetitions) was 10.338 seconds when using ACML compared with 18.184 seconds when using the managed implementation. The improvement of using an optimized linear algebra library is obvious.

Using Probability and Statistics APIs

The support for probability and statistics is one distinguishing features of the Math.NET Numerics library. The library contains a wide range of probability distributions, samplers, as well as functions for descriptive statistics.

Working with Probability Distributions

Math.NET Numerics implements a large number of probability distributions, both continuous (represented using the IContinuousDistribution interface) and discrete (the IDiscreteDistribution interface). There are about thirty different univariate distributions and six multivariate distributions.

Each distribution is represented by an object that has a constructor to create a distribution using the specified parameters. The distribution object implements one of the interfaces and provides properties to get the basic statistics such as the mean and standard deviation of the distribution. Both interfaces also provide a method Sample that generates a random sample from this distribution. This is a very important operation required by many algorithms such as Markov chain Monte Carlo (MCMC) methods.

The following listing creates a normal distribution and uses it to generate two random samples:

> open MathNet.Numerics.Distributions;;

> let normal = new Normal(1.0, 4.0);;
val normal : Normal = Normal(Mean = 1, StdDev = 4)

> normal.Sample();;
val it : float = 1.765599352

> normal.Sample();;
val it : float = 4.678777683

The objects representing probability distributions can be found in the MathNet.Numerics.Distributions namespace. When creating a Normal distribution, the constructor expects a mean and a standard deviation as arguments. After creating the distribution, the Sample method can be used to generate random samples. The next snippet shows how to obtain basic statistics about the distribution:

> normal.Density(5.0);;
val it : float = 0.06049268113

> normal.InverseCumulativeDistribution(0.95);;
val it : float = 7.579414508

> normal.Samples() |> Seq.take 10 |> Seq.toList;;
val it : float list =
    [ -0.6614167747; 2.742996449; 1.136350009; 
       2.081772036; 4.896975162; -3.450881607; 
       4.792252119; 1.316237836; -5.052077039; 1.461782835 ]

The first line shows the density value at 5.0 for normal. The second line shows how to get the right point value (7.579) given a cumulative sum (0.95) of a normal distribution. Finally, the Samples member can be used to generate an infinite sequence of random samples. The snippet uses Seq.take and Seq.toList to take first 10 samples and convert them to an F# list.

Using Descriptive Statistics

Another feature that the library provides is a set of extension methods for calculating descriptive statistics of data collections such as mean, variance, maximum, and minimum. The methods are implemented in a Statistics class in the MathNet.Numerics.Statistics namespace. They can be called either directly or using the dot-notation on any type that implements IEnumerable<float>:

> open MathNet.Numerics.Statistics;;
> let data = normal.Samples() |> Seq.take 100 |> Seq.toList;;
val data : float list = [ ... ]

> data.Mean();;
val it : float = 0.8555745637

> data.StandardDeviation();;
val it : float = 3.657119715

The snippet uses the normal distribution from the previous section to generate a random sequence of numbers using a normal distribution with mean 1.0 and a standard deviation 4.0. As expected, the statistics calculated from the data using the Mean and StandardDeviation extension methods are quite close to the properties of the distribution.

In addition, the library also provides a DescriptiveStatistics object that can be created from a collection and precomputes all the statistics at once:

> let stats = new DescriptiveStatistics(data);;
val stats : DescriptiveStatistics

> stats.Mean;;
val it : float = 0.8555745637

> stats.StandardDeviation;;
val it : float = 3.657119715

The listing uses the DescriptiveStatstics type to calculate the mean, standard deviation as well as numerous other statistics from the data. Since the snippet uses the same dataset, the results are the same as the values calculated using extension methods. The main benefit of using the DescriptiveStatstics type is that it processes the data only once.

Summary

This tutorial introduced various components of the Math.NET Numerics library. It discussed how to work with matrices and vectors using the .NET API provided by the library as well as using F# wrappers. The tutorial also demonstrated how to use linear algebra functionality of the library. The basic implementation is fully managed, but it is possible to use native libraries such as Intel MKL or AMD ACML.

Finally, the article looked at several topics from the area of probability and statistics. It demonstrated how to construct and sample a distribution and how to calculate basic descriptive statistics. Some of the topics that were not covered, but are supported by Math.NET Numerics, include numerical integration and interpolation.

Additional Resources

This article introduced the Math.NET Numerics library and explained how to use it from F#. Math.NET is just one of several libraries that can be used to write numerical computations. The following articles review several other options and also discuss visualization of data in F#:

• Tutorial: Using the F# PowerPack and the F# MathProvider

• Using Microsoft Sho in F#

• Overview: Charting Libraries for F#

• Overview: Getting Started with the FSharpChart Library

As a functional programming language, F# differs in many ways from languages like R, Matlab, or Python. The following articles discuss the key functional concepts and their effect when implementing parallel computations:

• Overview: Introducing F# and Functional Programming

• Overview: Parallel and Asynchronous Programming

The following documents and online resources provide more information about Math.NET Numerics and the AMD ACML library that was used in this tutorial:

• Math.NET Numerics is the project homepage containing downloads, documentation, and other resources about the library.

• AMD Core Math Library (ACML) contains more information about the AMD ACML library as well as downloads for various platforms.

To download the code snippets shown in this article, go to https://code.msdn.microsoft.com/Chapter-4-Numerical-3df3edee

See Also

This article is based on Real World Functional Programming: With Examples in F# and C#. Book chapters related to the content of this article are:

• Book Chapter 10: “Efficiency of data structures” explains how to write efficient programs in F#. The chapter covers advanced functional techniques and explains how and when to use arrays and functional lists.

• Book Chapter 12: “Sequence expressions and alternative workflows” contains detailed information on processing in-memory data (such as lists and seq<'T> values) using higher-order functions and F# sequence expressions.

• Book Chapter 13: “Asynchronous and data-driven programming” explains how asynchronous workflows work and uses them to write an interactive script that downloads a large dataset from the Internet.

Previous article: Using Microsoft Sho in F#

Next article: Tutorial: Using the F# PowerPack and the F# MathProvider