What's New in Math.NET Numerics 2.6 | Christoph Rüegg

What's New in Math.NET Numerics 2.6

Math.NET Numerics v2.6, released in July 2013, is focused on filling some gaps around the very basic numerical problems of fitting a curve to data and finding solutions of nonlinear equations. As usual you'll find a full listing of all changes in the release notes. However, I'd like to take the chance to highlight some important changes, show some code samples and explain the reasoning behind the changes.

A lot of high quality code contributions made this release possible. Just like last release, I've tried to attribute them directly in the release notes. Thanks again!

Please let me know if these "What's New" articles are useful in this format and whether I should continue to put the together for future releases. See also what's new in the previous version 2.5.

Linear Curve Fitting

Fitting a linear-parametric curve to a set of samples such that the squared errors are minimal has always been possible with the linear algebra toolkit, but it was somewhat complicated to do and required understanding of the algorithm. See Linear Regression with Math.NET Numerics for an introduction and some examples.

Note: if you need to have the curve go exactly through all your data points, use our Interpolation routines instead.

We now finally provide a shortcut with a few common functions to fit to data, but also a method to fit a linear combination of arbitrary functions. For fitting a simple line it uses an efficient direct algorithm:

1: 
2: 
3: 
4: 
5: 
6: 
7: 
8: 
var x = new [] { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 };
var y = new [] { 4.986, 2.347, 2.061, -2.995, 2.352, -5.782 }

C#: var p = Fit.Line(x, y);
    var offset = p[0];  // = 7.01013
    var slope = p[1];   // = -2.08551

F#: let offset, slope = Fit.line x y

Otherwise it usually applies an ordinary least squares regression to find the best parameters using a thin QR decomposition (leveraging a native provider like Intel MKL if enabled). This also works with arbitrary functions, like sine and cosine:

1: 
2: 
3: 
4: 
F#: let p = (x, y) ||> Fit.linear [(fun _ -> 1.0); (Math.Sin); (Math.Cos)]
C#: var p = Fit.LinearCombination(x, y, z => 1.0, Math.Sin, Math.Cos);

// p = [ -0.287, 4.02, -1.46 ], hence f: x -> -0.287 + 4.02*sin(x) - 1.46*cos(x)

The intention is to add more special cases for common curves like the logistic function in the future. Like the line they may have more appropriate direct implementations. For now there is one other special case, for fitting to a polynomial. It returns the best parameters, in ascending order (coefficient for power k has index k) compatible to the Evaluate.Polynomial routine:

1: 
2: 
3: 
C#: var coeff = Fit.Polynomial(x, y, 2); // order 2
    Evaluate.Polynomial(1.2, coeff); // ...
F#: let coeff = Fit.polynomial 2 x y

In practice your x values are not always just real numbers. Maybe you need multi-dimensional fitting where the x values are actually arrays, or even full data structures. For such cases we provide a version that is generic in x and where you can provide a list of functions that accept such x directly without the need to convert to an intermediate double vector first:

1: 
2: 
3: 
C#: var p = Fit.LinearMultiDim(xarrays, y, f1, f2, f3, ...);
    var p = Fit.LinearGeneric(xstructs, y, f1, f2, f3, ...);
F#: let p = Fit.linear [f1; f2; f3; ...] xgeneric y

Often after evaluating the best fitting linear parameters you'd actually want to evaluate the function with those parameters. For this scenario we provide a shortcut as well: For each of these methods there is also a version with a "Func" suffix ("F" in F#) which, instead of the parameters, returns the composed function:

1: 
2: 
3: 
4: 
5: 
F#: let f = Fit.lineF x y
    [1.0..0.1..2.0] |> List.map f

C#: var f = Fit.LinearCombinationFunc(x, y, z => z*z, Math.Sin, SpecialFunctions.Gamma);
    Enumerable.Range(0,11).Select(x => f(x/10.0))

Root Finding

We now provide basic root finding algorithms. A root of a function x -> f(x) is a solution of the equation f(x)=0. Root-finding algorithms can thus help finding numerical real solutions of arbitrary equations, provided f is reasonably well-behaved and we already have an idea about an interval [a,b] where we expect a root. As usual, there is a facade class FindRoots for simple scenarios:

The routines usually expect a lower and upper boundary as parameters, and then optionally the accuracy we try to achieve and the maximum number of iterations.

1: 
2: 
3: 
4: 
5: 
6: 
C#: FindRoots.OfFunction(x => x*x - 4, -5, 5) // -2.00000000046908
C#: FindRoots.OfFunction(x => x*x - 4, -5, 5, accuracy: 1e-14) // -2 (exact)
C#: FindRoots.OfFunctionDerivative(x => x*x - 4, x => 2*x, -5, 5) // -2 (exact)

F#: FindRoots.ofFunction -5.0 5.0 (fun x -> x*x - 4.0)
F#: FindRoots.ofFunctionDerivative -5.0 5.0 (fun x -> x*x - 4.0) (fun x -> 2.0*x)

A NonConvergenceException is thrown if no root can be found by the algorithm.

In practice you'd often want to use a specific well-known algorithm. You'll find them in the RootFinding namespace. Each of these algorithm provides a FindRoot method with similar arguments as those above. However, the algorithms may sometimes fail to find a root or the function may not actually have a root within the provided interval. Failing to find a root is thus not exactly exceptional. That's why the algorithms also provide an exception-free TryFindRoot code path with the common Try-pattern as in TryParse.

Bisection

A simple and robust yet rather slow algorithm, implemented in the Bisection class.

Example: Find the real roots of the cubic polynomial 2x^3 + 4x^2 - 50x + 6:

Root Finding Equation 2

1: 
2: 
3: 
4: 
5: 
6: 
Func<double, double> f = x => Evaluate.Polynomial(x, 6, -50, 4, 2);
Bisection.FindRoot(f, -6.5, -5.5, 1e-8, 100);  // -6.14665621970684
Bisection.FindRoot(f, -0.5, 0.5, 1e-8, 100);   // 0.121247371972135
Bisection.FindRoot(f, 3.5, 4.5, 1e-8, 100);    // 4.02540884774855

F#: f |> FindRoots.bisection 100 1e-8 3.5 4.5  // Some(4.0254..)

Note that the F# function returns a float option. Instead of throwing an exception it will simply return None if it fails.

Brent's Method

We use Brent's method as default algorithm, implemented in the Brent class. Brent's method is faster than bisection, but falls back to something close to bisection if the faster approaches (essentially the secant method and inverse quadratic interpolation) fails and is therefore almost as reliable.

The same example as above, but using Brent's method:

1: 
2: 
3: 
4: 
5: 
6: 
Func<double, double> f = x => Evaluate.Polynomial(x, 6, -50, 4, 2);
Brent.FindRoot(f, -6.5, -5.5, 1e-8, 100);  // -6.14665621970684
Brent.FindRoot(f, -0.5, 0.5, 1e-8, 100);   // 0.121247371972135
Brent.FindRoot(f, 3.5, 4.5, 1e-8, 100);    // 4.02540884774855

F#: f |> FindRoots.brent 100 1e-8 3.5 4.5  // Some(4.0254..)

Note that there are better algorithms for finding all roots of a polynomial. We plan to add specific polynomial root finding algorithms later on.

Newton-Raphson

The Newton-Raphson method leverages the function's first derivative to converge much faster, but can also fail completely. The pure Newton-Raphson algorithm is implemented in the NewtonRaphson class. However, we also provide a modified algorithm that tries to recover (instead of just failing) when overshooting, converging too slowly or even when loosing bracketing in the presence of a pole. This modified algorithms is available in the RobustNewtonRaphson class.

Example: Assume we want to find solutions of x+1/(x-2) == -2, hence x -> f(x) = 1/(x-2)+x+2 with a pole at x==2:

Root Finding Equation 3

1: 
2: 
3: 
4: 
5: 
6: 
7: 
8: 
9: 
Func<double, double> f = x => 1/(x - 2) + x + 2;
Func<double, double> df = x => -1/(x*x - 4*x + 4) + 1;
RobustNewtonRaphson.FindRoot(f, df, -2, -1, 1e-14, 100, 20); // -1.73205080756888
RobustNewtonRaphson.FindRoot(f, df, 1, 1.99, 1e-14, 100, 20); // 1.73205080756888
RobustNewtonRaphson.FindRoot(f, df, -1.5, 1.99, 1e-14, 100, 20); // 1.73205080756888
RobustNewtonRaphson.FindRoot(f, df, 1, 6, 1e-14, 100, 20); // 1.73205080756888

F#: FindRoots.newtonRaphsonRobust 100 20 1e-14 1.0 6.0 f df
F#: (f, df) ||> FindRoots.newtonRaphsonRobust 100 20 1e-14 1.0 6.0

Broyden's Method

The quasi-newton method by Broyden, implemented in the Broyden class, may help you to find roots in multi-dimensional problems.

Linear Algebra

As usual there have been quite a few improvements around linear algebra, see the release notes for the complete list. If you've enabled our Intel MKL native linear algebra provider, then eigenvalue decompositions should be much faster now. Matrices now also support the new F# 3.1 array slicing syntax.

Note that we're phasing out the MathNet.Numerics.IO library and namespace and plan to drop it entirely in v3. We've already replaced it with two new separate NuGet packages and obsoleted all members of the old library. The new approach with separate libraries makes it possible to introduce specific dependencies e.g. to read and write Excel files, without forcing these dependencies on all of Math.NET Numerics. We recommend to switch over to the new packages as soon as possible.

Statistics

We've had a Pearson correlation coefficient routine for a while, but no Covariance routine. In addition to a new Spearman ranked correlation routine, this release finally also adds sample and population Covariance functions for arrays and IEnumerables.

1: 
ArrayStatistics.Covariance(new[] {1.2, 1.3, 2.4}, new[] {2.2, 2.3, -4.5})
Multiple items
module List

from Microsoft.FSharp.Collections

--------------------
type List<'T> =
  | ( [] )
  | ( :: ) of Head: 'T * Tail: 'T list
  interface IEnumerable
  interface IEnumerable<'T>
  member GetSlice : startIndex:int option * endIndex:int option -> 'T list
  member Head : 'T
  member IsEmpty : bool
  member Item : index:int -> 'T with get
  member Length : int
  member Tail : 'T list
  static member Cons : head:'T * tail:'T list -> 'T list
  static member Empty : 'T list

Full name: Microsoft.FSharp.Collections.List<_>
val map : mapping:('T -> 'U) -> list:'T list -> 'U list

Full name: Microsoft.FSharp.Collections.List.map
Multiple items
val double : value:'T -> double (requires member op_Explicit)

Full name: Microsoft.FSharp.Core.ExtraTopLevelOperators.double

--------------------
type double = System.Double

Full name: Microsoft.FSharp.Core.double