
Hello,
I am currently running numerical simulations and have been using F# for the task. F# has a flexible syntax that allows for expressive mathematical notation as well as having decent performance (with better performance possible via dropping into C++/CLI
as needed). I like the clean OO API being used in Numerics, and I think that definitely is the way to go for the underlying numerical library. However, I was wondering if there would be any interest in creating a more "scriptlike" environment for
the F# bindings so as to allow clean mathematical notation? A good example of this approach is the SciPy/NumPy/Matplotlib stack, where the core algorithms are in an OO hierarchy but they also expose a scripting environment called pylab which lets you write
code in a quick and clean style similar to MATLAB/Octave.
There are really two parts to such an environment: more overloaded operators on the vector/matrix types and a library of functions which follows MATLABlike naming conventions and allow easy generation of numerical code without getting into the depths of
the OO hierarchy. For example, it'd be nice to be able to write a=zeros(3,3) and have a 3x3 matrix of zeros, without figuring out what storage type you want (as mentioned in a previous thread, 95% of the time you probably want Double.DenseMatrix or Complex.DenseMatrix).
To that end I wanted to suggest reforming the F# bindings currently in the project to take advantage of the F# language's flexibility in creating a scripting environment.
Here is what some example code would look like:
let a = zeros(3,3) // 3x3 zero matrix
let b = ones(3,3) // 2x2 matrix of ones
let c = b.^3 // c=b*b*b
let d = b.^(1) // d=inverse(b) (should probably be pseudoinv too?)
let e = b.*b // e is all ones, elementwise multiplication
let f = b*b // f is all threes, matrix multiplication
let g = b*b*2.0+b./2.0 // composition works fine
The choice of operators can exactly copy MATLAB, since periods are usable as operators in F#. Here is a toy example of what a Kalman Filter would look like using such a binding:
let KalmanFilter (model:KFModel) (estimate:KFEstimate) (z:Matrix)=
let (x,P)=(estimate.x,estimate.P)
let (H,R,Phi,Q)=(model.H,model.R,model.Phi,model.Q)
// First propagate to the current time
let x = Phi*x
let P = Phi*P*Phi.T+Q
// Now update using measurement
let K = P*H.T*(H*P*H.T+R).^(1)
let res = zH*x
let x = x+K*res
let P = P  K*H*P
The code is clean and readable by a mathematician/engineer without parsing parentheses. Of course, talk is cheap, so here is a partial implementation of a Matrix<T> wrapper type in F# which would allow the above code to run. I believe there is no easy
way to do this without wrapping the C# class in a F# type unless the C# implementation of Matrix<T> is modified, as type extensions do not allow operator overloading and static type parameters aren't flexible enough. Of course, I'd love to hear a better
way! Hopefully this is the right place for this, I wasn't sure what the codeplex conventions were for code pasting:
type Matrix<'T when 'T : (new:unit > 'T) and 'T:struct and 'T :>System.ValueType and 'T :> System.IEquatable<'T> and 'T :> System.IFormattable> (inner:Generic.Matrix<'T>) =
member this.inner with get() = inner
member this.Item
with get(n,m)=inner.[n,m]
and set(n,m) v=inner.[n,m]<v
override this.ToString()=
inner.ToString()
// Redefine transpose for brevity
member this.T with get()=new Matrix<'T>(inner.Transpose())
// Hermitian transpose
member this.H with get()=new Matrix<'T>(inner.ConjugateTranspose())
// Elementwise operations (matching MATLAB syntax)
static member (.*) (a:Matrix<'T>,b:Matrix<'T>) = new Matrix<'T>(a.inner.PointwiseMultiply(b.inner))
static member (./) (a:Matrix<'T>,b:Matrix<'T>) = new Matrix<'T>(a.inner.PointwiseDivide(b.inner))
// Elementwise work on scalars (following MATLAB convention)
static member (.*) (a:Matrix<'T>,b:'T) = new Matrix<'T>(a.inner.Multiply(b))
static member (./) (a:Matrix<'T>,b:'T) = new Matrix<'T>(a.inner.Divide(b))
/// Matrix scalar exponentiation for integer powers > 0, inverses
static member (.^) (a:Matrix<'T>,b:int) =
if b= 1 then
new Matrix<'T>(a.inner.Inverse())
elif b<1 then
raise (new NotImplementedException())
else
let mutable out = a.inner.Clone()
for i in 0..b2 do
out<out*a.inner
new Matrix<'T>(out)
/// Matrixvalued exponentiation (i.e. e^(matrix))
static member (.^) (a:double,b:Matrix<'T>) =
raise (new NotImplementedException()) // todo
/// Matrix "right" division (i.e. conventional division)
static member (/) (a:Matrix<'T>,b:Matrix<'T>) =
if a.inner.ColumnCount=a.inner.RowCount && b.inner.ColumnCount=b.inner.RowCount && b.inner.Rank()=b.inner.ColumnCount then
new Matrix<'T>(a.inner*b.inner.Inverse())
else
raise (new NotImplementedException()) // todo: Should use pseudoinverse
/// Matrix "left" division (i.e. mathemtically this performs a\b, but "\" is not a valid F# operator)
static member ldiv (a:Matrix<'T>,b:Matrix<'T>) =
if a.inner.ColumnCount=a.inner.RowCount && b.inner.ColumnCount=b.inner.RowCount && a.inner.Rank()=a.inner.ColumnCount then
new Matrix<'T>(a.inner.Inverse()*b.inner)
else
raise (new NotImplementedException()) // todo: Should use pseudoinverse
// Conventional operations. These are already overloaded in C# for the inner class
static member (*) (a:Matrix<'T>,b:Matrix<'T>) = new Matrix<'T>(a.inner.Multiply(b.inner))
static member (*) (a:Matrix<'T>,b:'T) = new Matrix<'T>(a.inner.Multiply(b))
static member (*) (a:'T,b:Matrix<'T>) = new Matrix<'T>(b.inner.Multiply(a))
static member () (a:Matrix<'T>,b:Matrix<'T>) = new Matrix<'T>(a.inner.Subtract(b.inner))
static member (~) (a:Matrix<'T>) = new Matrix<'T>(a.inner.Negate())
static member (+) (a:Matrix<'T>,b:Matrix<'T>) = new Matrix<'T>(a.inner.Add(b.inner))
static member (~+) (a:Matrix<'T>) = a
We'd also need the supporting library functions we used:
let zeros(n,m)=
new Matrix(new Double.DenseMatrix(n,m))
let ones(n,m)=
new Matrix(new Double.DenseMatrix(n,m,1.0))
let linspace(a,b,length)=
let mutable out = zeros(1,length)
for i in 0..length1 do
out.[0,i]<a+double(i)*(ba)/double(length1)
out
If this is something that interests anyone, I could start writing a more complete wrapper that wraps all the exposed members of Matrix<T>. Also, the current F# bindings (i.e. ofList, etc.) would need to be modified to support the F# wrapper class,
as they are quite useful and idiomatic F#. Let me know what everyone thinks, and thanks for actually reading my giant wall of text/code! :)
P.S. while typing up this example I think I may have run into a bug. Apparently the code:
let g = MathNet.Numerics.LinearAlgebra.Double.DenseMatrix(2,2,3.0).Divide(2.0)
returns a zero matrix. Is that behavior intended?



I would like to second this idea. I'm not sure about implementation or anything, but F# seems ideally suited for numerics, so an F# interface would be great. I can contribute cycles if needed.


Mar 31, 2011 at 3:30 PM
Edited Mar 31, 2011 at 3:47 PM

>let g = MathNet.Numerics.LinearAlgebra.Double.DenseMatrix(2,2,3.0).Divide(2.0) >returns a zero matrix. Is that behavior intended?
I'll look in to it. Thanks for pointing it out.
As for your suggestion, I like the idea but don't know enough about F# to really comment. Jurgen, what do you think?



This discussion has been copied to a work item. Click
here to go to the work item and continue the discussion.



>let g = MathNet.Numerics.LinearAlgebra.Double.DenseMatrix(2,2,3.0).Divide(2.0)
>returns a zero matrix. Is that behavior intended?
I wasn't able to duplicate the problem. Please try the most recent change set.



Yeah my bad, all is well if I compile from source (was using beta 1).



I've had some time to improve the F# interface I posted earlier:
https://gist.github.com/922192
There are two types defined here:
 FSMatrix is a wrapper around Matrix<T> that has the usual overloaded operators
 NumericsPlot is a class which exposes a simple plotting interface. You can pass it X and Y coordinates or just Y coordinates (X will be 1,2,3...). It accepts a list of points, array of points, Matrix<T> of points, or FSMatrix of points. Only works for
integers/floats. When plotting matrices it plots each column as a separate data set. Note that plots are run in their own threads with their own event loops, so your program will not halt until the windows are closed (therefore if you want to just plot some
stuff, you can exit after calling plot() and the program will stay open).
Here is an example of plotting a matrix/simple list:
open System
open FSNumerics
open MathNet.Numerics.LinearAlgebra
// Plot a list of points
let p = new NumericsPlot()
p.plot([1;2;3;4;5;4;3;2;1])
// Do an experiment, fill a matrix, print the results
let numRuns = 5
let numSamples = 50000
let data = zeros(numRuns,numSamples)
for i in 0..data.inner.RowCount1 do
for j in 0..data.inner.ColumnCount1 do
data.[i,j]<sin(float(j)/10000.0float(i*2))
p.plot(data.T) // Plots each column in the matrix in a separate plot, total of 5 plots
The resulting plots look like this:
http://i.imgur.com/hyeHp.png



Hi kyon,
I like the succinct syntax of this new type. It's too bad we can use F# type augmentations to augment the existing types with the operators you've defined. I think that would be a neater solution rather than defining a whole new type. I'm not entirely sure
we should add yet another specific F# type to the library; let me get in touch with some other people from the F# community and have them weigh in ...
(Sorry for the massive delay)



Hi!
Jurgen asked me for some thoughts about the suggestion, so I played a bit with using Math.NET from F#. The wrapper definitely makes some common tasks nicer and there isn't really a way to get exactly the same extensions without using a wrapper. On the other
hand, wrappers are (in principle) a bit painful. If you want to wrap everything, then it is a lot of code that must be maintained and if it is not complete, then people will need to explicitly wrap/unwrap the matrix themselves to get all the features (for
example, the extension methods for Svd, Evd, GramSchmidt, etc. would have to be duplicated too).
You can get some of the features just by adding extension members to the type in the F# library. For example, here is a definition of the T and H members:
type MathNet.Numerics.LinearAlgebra.Generic.
Matrix<'T when 'T : struct and 'T : (new : unit > 'T) and 'T :> System.ValueType and
'T :> System.IEquatable<'T> and 'T :> System.IFormattable> with
// Redefine transpose for brevity
member this.T = this.Transpose()
// Hermitian transpose
member this.H = this.ConjugateTranspose()
The additional operators are more tricky  they can be defined using "let" outside of the type, but then they cannot be overloaded (which means you'd need different operator for pointwise multiplication and multiplication by a scalar). One way to implement
those is to use static member constraints, which essentially just says  if the type has "PointwiseMultiply" method, then it can use this operator:
let inline (.*) (a:^T) (b:^T) =
(^T : (member PointwiseMultiply : ^T > ^R) (a, b))
let inline (./) (a:^T) (b:^T) =
(^T : (member PointwiseDivide : ^T > ^R) (a, b))
let inline (.*) (a:^T) (b:^S) =
(^T : (member Multiply : ^S > ^R) (a, b))
let inline (./) (a:^T) (b:^S) =
(^T : (member Divide : ^S > ^R) (a, b))
The nice thing about this trick is that it will work for both matrices and vectors. You could even use the same operator for different types of operations (e.g. scalar and pointwise multiplication), but then the Matrix type would have to provide these two
operations as some overloaded method with the same name (it could be something like "FSharpOverloadedMultiply", but that would look ugly from C#). It is not as nice as what the wrapper can provide, but it looks decent:
let zeros n m = new Double.DenseMatrix(n,m)
let ones n m = new Double.DenseMatrix(n,m,1.0)
let v = vector [ 1.0; 2.0 ]
v .* (v .* 2.0)
let m = matrix [ [ 1.0; 2.0]; [ 3.0; 4.0 ] ]
m.H .* ((ones 2 2) .* 2.0)
Regarding charting  there is an
F# sample called FSharpChart that implements a pretty complete F# wrapper for the .NET WinForms charting libraries. The code is available as a sample with permissive license, so it should be fine to include it in the F# library for Math.NET Numerics (or
in some other part of Math.NET). That would be actually quite nice, because currently it is "just" a sample.
I don't have any clear suggestion on what is better design  adding F# specific functionality using extensions and additional operators doesn't give you all the features that a wrapper can give, but wrappers have quite a few issues. I guess it is also perfectly
fine to have the wrapper in some additional script file that people can reference if they want, but maybe not in the standard library.
Hope this helps!
Tomas

