.net - F# 的简单包装器，用于执行矩阵运算

Question

这是一篇比较长的文章。 F# 有一个 matrix和 vector现在键入(在 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"

我也和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

时序(双核 2.0GH，2GB 内存，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还有similar results .

Answer 1

我认为 dnmat2mat 和 randmat 都可以通过使用 Matrix.init 来简化:

let dnmat2mat (dnmat : DenseMatrix) =
  Matrix.init (dnmat.Rows) (dnmat.Columns) (fun i j -> dnmat.[i,j])

let randmat n m =
  let r = System.Random()
  Matrix.init n m (fun _ _ -> r.NextDouble())