A simple wrapper for F # to perform matrix operations

Question

This is a relatively long post. F # has a matrix and vector type (not in the kernel in PowerPack). Fine! Even the numerical computational power of Python belongs to the third part.

But the functions provided there are limited by matrix and vector arithmetic, therefore, for performing inversions, decompositions, etc. we still need to use another library. Now I use the latest dnAnalytics , which is integrated into the Math.Net project. But the Math.Net project has not had any updates to the public for more than a year, I don’t know if they have a continuation plan.

I made the following shell, this shell uses Matlab-like functions for simple linear algebra. Since I am new to F # and FP, could you give some tips on improving the shell and code? Thanks!

#r @"D:\WORK\tools\dnAnalytics_windows_x86\bin\dnAnalytics.dll" #r @"FSharp.PowerPack.dll" open dnAnalytics.LinearAlgebra open Microsoft.FSharp.Math open dnAnalytics.LinearAlgebra.Decomposition // F# matrix -> ndAnalytics DenseMatrix let mat2dnmat (mat:matrix) = let m = new DenseMatrix(mat.ToArray2D()) m // ndAnalytics DenseMatrix -> F# matrix let dnmat2mat (dnmat:DenseMatrix) = let n = dnmat.Rows let m = dnmat.Columns let mat = Matrix.create nm 0. for i=0 to n-1 do for j=0 to m-1 do mat.[i,j] <- dnmat.Item(i,j) mat // random matrix let randmat nm = let r = new System.Random() let ranlist m = [ for i in 1..m do yield r.NextDouble() ] matrix ([1..n] |> List.map (fun x-> ranlist m)) // is square matrix let issqr (m:matrix) = let n, m = m.Dimensions n = m // is postive definite let ispd m = if not (issqr m) then false else let m1 = mat2dnmat m let qrsolver = dnAnalytics.LinearAlgebra.Decomposition.Cholesky(m1) qrsolver.IsPositiveDefinite() // determinant let det m = let m1 = mat2dnmat m let lusolver = dnAnalytics.LinearAlgebra.Decomposition.LU(m1) lusolver.Determinant () // is full rank let isfull m = let m1 = mat2dnmat m let qrsolver = dnAnalytics.LinearAlgebra.Decomposition.GramSchmidt(m1) qrsolver.IsFullRank() // rank let rank m = let m1 = mat2dnmat m let svdsolver = dnAnalytics.LinearAlgebra.Decomposition.Svd(m1, false) svdsolver.Rank() // inversion by lu let inv m = let m1 = mat2dnmat m let lusolver = dnAnalytics.LinearAlgebra.Decomposition.LU(m1) lusolver.Inverse() // lu let lu m = let m1 = mat2dnmat m let lusolver = dnAnalytics.LinearAlgebra.Decomposition.LU(m1) let l = dnmat2mat (DenseMatrix (lusolver.LowerFactor ())) let u = dnmat2mat (DenseMatrix (lusolver.UpperFactor ())) (l,u) // qr let qr m = let m1 = mat2dnmat m let qrsolver = dnAnalytics.LinearAlgebra.Decomposition.GramSchmidt(m1) let q = dnmat2mat (DenseMatrix (qrsolver.Q())) let r = dnmat2mat (DenseMatrix (qrsolver.R())) (q, r) // svd let svd m = let m1 = mat2dnmat m let svdsolver = dnAnalytics.LinearAlgebra.Decomposition.Svd(m1, true) let u = dnmat2mat (DenseMatrix (svdsolver.U())) let w = dnmat2mat (DenseMatrix (svdsolver.W())) let vt = dnmat2mat (DenseMatrix (svdsolver.VT())) (u, w, vt.Transpose) 

and test code

 (* todo: read matrix market format ref: http://math.nist.gov/MatrixMarket/formats.html *) let readmat (filename:string) = System.IO.File.ReadAllLines(filename) |> Array.map (fun x-> (x |> String.split [' '] |> List.toArray |> Array.map float)) |> matrix let timeit f str= let watch = new System.Diagnostics.Stopwatch() watch.Start() let res = f() watch.Stop() printfn "%s Needed %f ms" str watch.Elapsed.TotalMilliseconds res let test() = let testlu() = for i=1 to 10 do let a,b = lu (randmat 1000 1000) () () let testsvd() = for i=1 to 10 do let u,w,v = svd (randmat 300 300) () () let testdet() = for i=1 to 10 do let d = det (randmat 650 650) () () timeit testlu "lu" timeit testsvd "svd" timeit testdet "det" 

I also compared with matlab

 t = cputime; for i=1:10, [l,u] = lu(rand(1000,1000)); end; e = cputime-t t = cputime; for i=1:10, [u,w,vt] = svd(rand(300,300)); end; e = cputime-t t = cputime; for i=1:10, d = det(rand(650,650)); end; e = cputime-t 

Timings (Duo Core 2.0GH, 2 GB memory, Matlab 2009a)

 f#: lu Needed 8875.941700 ms svd Needed 14469.102900 ms det Needed 2820.304600 ms matlab: lu 3.7752 svd 5.7408 det 1.2636 

matlab is about twice as fast. This is reasonable since the native R also has similar results .

Answer 1

I think that both dnmat2mat and randmat could be simplified using Matrix.init :

 let dnmat2mat (dnmat : DenseMatrix) = Matrix.init (dnmat.Rows) (dnmat.Columns) (fun ij -> dnmat.[i,j]) let randmat nm = let r = System.Random() Matrix.init nm (fun _ _ -> r.NextDouble())