Простая оболочка для F# для выполнения матричных операций

Это относительно длинный пост. F# имеет матрица и векторный тип (в PowerPack, а не в Core) теперь. Это здорово! Даже способность Python к численным вычислениям взята из третьей части.

Но функции, представленные там, ограничены матричной и векторной арифметикой, поэтому для инверсии, разложения и т. д. нам все равно нужно использовать другую библиотеку. Сейчас я использую последнюю версию dnAnalytics, которая объединяется с проектом Math.Net. Но у проекта Math.Net нет обновлений для публики уже больше года, я не знаю, есть ли у них планы на продолжение.

Я сделал следующую оболочку, эта оболочка использует функции, подобные Matlab, для выполнения простой линейной алгебры. Поскольку я новичок в F # и FP, не могли бы вы дать несколько советов по улучшению оболочки и кода? Спасибо!

#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 n m 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 n m =
    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)

и тестовый код

(* 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"

я тоже сравнивал с матлабом

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

Тайминги (Duo Core 2.0GH, 2 Гб памяти, 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 примерно в два раза быстрее. Это разумно, так как собственный R также имеет похожие результаты.