Un wrapper simple pour F # pour faire des opérations matricielles

Question

Il s'agit d'un article relativement long. F # a un type matrix et vector (dans PowerPack pas dans le Core) maintenant. C'est bien! Même la capacité de calcul numérique de Python provient de la troisième partie.

Mais les fonctions fournies ici sont limitées à l'arithmétique matricielle et vectorielle, donc pour faire l'inversion, les décompositions, etc., nous avons encore besoin d'utiliser une autre bibliothèque. J'utilise maintenant la dernière dnAnalytics, qui fusionne dans le projet Math.Net. Mais le projet Math.Net n'a pas de mises à jour pour le public depuis plus d'une année maintenant, je ne sais pas si elles ont un plan pour continuer.

J'ai fait le wrapper suivant, ce wrapper utilise des fonctions de type Matlab pour faire de l'algèbre linéaire simple. Comme je suis nouveau à F # et FP, pourriez-vous donner quelques conseils pour améliorer l'emballage et le code? Merci!

#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) 

et code de test

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

J'ai aussi comparé à 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 

Les timings (Core Duo 2.0GH, 2 Go de mémoire, 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 est environ deux fois plus rapide. Ceci est raisonnable, car R natif a également similar results.

Answer 1

Je pense que les deux dnmat2mat et randmat pourrait être simplifié en utilisant 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())