Vector, Matrix and Linear Algebra in F#

The following script snippets present the vector and matrix functionalities in Microsoft.Fsharp.Numerics namespace from Fsharp.PowerPack assemble.
Also, they present the Linear Algebra functionalities from FlyingFrog namespace from FSharpForNumeric assemble.

//F# documentation at
//http://msdn.microsoft.com/library/dd233154(VS.100).aspx
//Fsharp.PowerPack documentation at
//http://research.microsoft.com/en-us/um/cambridge/projects/fsharp/manual/namespaces.html
//Flying frogs numerics documentation at
//http://www.ffconsultancy.com/products/fsharp_for_numerics/docs/namespaces.html


#I @"C:\Program Files (x86)\FSharpPowerPack-2.0.0.0\bin"
#r "Fsharp.PowerPack.dll"
#r "FSharpForNumerics.dll"
#r "FSharpForVisualization.dll"

open Microsoft.FSharp.Math
open FlyingFrog
open FlyingFrog.Graphics
open FlyingFrog.Graphics.Typeset;;

//vector (float)
let v1 = vector [1.0;1.0;1.0]
let v2 = vector [2.0;2.0;2.0]
let v3 = rowvec [3.0;3.0;3.0];;val v1 : vector = vector [|1.0; 1.0; 1.0|]
val v2 : vector = vector [|2.0; 2.0; 2.0|]
val v3 : rowvec = rowvec [|3.0; 3.0; 3.0|]

//Simple operations
2.0*v1;;
v1+v2;;
v1-v2;;
v1.*v2;; //per element multiplication
let m1 = v2*v3;;
val it : Vector = vector [|2.0; 2.0; 2.0|]
val it : Vector = vector [|3.0; 3.0; 3.0|]
val it : Vector = vector [|-1.0; -1.0; -1.0|]
val it : Vector = vector [|2.0; 2.0; 2.0|]
val m1 : Matrix = matrix [[6.0; 6.0; 6.0]
[6.0; 6.0; 6.0]
[6.0; 6.0; 6.0]]
//Transpose
v1.Transpose;;val it : RowVector = rowvec [|1.0; 1.0; 1.0|]

//Vector Generic
type intVec = Vector
let intVec ll = Vector.Generic.ofSeq ll
let A = intVec [1;2;3;4]
let B = intVec [1;1;1;1]
A+B;;
A*10;;
Vector.Generic.sum A;;
type intVec = Vector
val intVec : seq<'a> -> Vector<'a>
val A : Vector = vector [|1; 2; 3; 4|]
val B : Vector = vector [|1; 1; 1; 1|]

>val it : Vector = vector [|2; 3; 4; 5|]
> val it : Vector = vector [|10; 20; 30; 40|]
> val it : int = 10

//matrix (float)
let A = matrix [ [ 1.0; 7.0; 2.0 ];
[ 1.0; 3.0; 1.0 ];
[ 2.0; 9.0; 1.0 ]; ]
let B = matrix [ [ 10.0; 70.0; 20.0 ];
[ 10.0; 30.0; 10.0 ];
[ 20.0; 90.0; 10.0 ]; ];;
val A : matrix = matrix [[1.0; 7.0; 2.0]
[1.0; 3.0; 1.0]
[2.0; 9.0; 1.0]]
val B : matrix = matrix [[10.0; 70.0; 20.0]
[10.0; 30.0; 10.0]
[20.0; 90.0; 10.0]]
//print
let f (i:Matrix<_>)= printfn "%A" (i.ToArray2D())
let g (i:Vector<_>)= printfn "%A" i;;
val f : Matrix<'a> -> unit
val g : Vector<'a> -> unit

// sum, subtract
f (A+B)> [[11.0; 77.0; 22.0]
[11.0; 33.0; 11.0]
[22.0; 99.0; 11.0]]
> [[-9.0; -63.0; -18.0]
[-9.0; -27.0; -9.0]
[-18.0; -81.0; -9.0]]

// matrix product
f (A*B);;
[[120.0; 460.0; 110.0]
[60.0; 250.0; 60.0]
[130.0; 500.0; 140.0]]

// element-wise product
f (A.*B);;
[[10.0; 490.0; 40.0]
[10.0; 90.0; 10.0]
[40.0; 810.0; 10.0]]
// scalar product
f (A * 2.0);;
[[2.0; 14.0; 4.0]
[2.0; 6.0; 2.0]
[4.0; 18.0; 2.0]]
f (2.0 * A);;
[[2.0; 14.0; 4.0]
[2.0; 6.0; 2.0]
[4.0; 18.0; 2.0]]
// negation of a matrix
f (-A);;
[[-1.0; -7.0; -2.0]
[-1.0; -3.0; -1.0]
[-2.0; -9.0; -1.0]]

// matrix-vector product
let b = vector [5.;8.;9.];;
g (A*b);;
val b : vector = vector [|5.0; 8.0; 9.0|]
> vector [|79.0; 38.0; 91.0|]

// Dimension
let dim = A.Dimensions;;
val dim : int * int = (3, 3)

//Transpose
Let A' = A.Transpose;;
val A' : Matrix = matrix [[1.0; 1.0; 2.0]
[7.0; 3.0; 9.0]
[2.0; 1.0; 1.0]]
//Number of rows,columns
let nrow = A.NumRows
let ncol = A.NumColslet;;
val nrow : int = 3
val ncol : int = 3

//Copy
let Anew = A.Copy();;
val Anew : Matrix = matrix [[1.0; 7.0; 2.0]
[1.0; 3.0; 1.0]
[2.0; 9.0; 1.0]]
// convert to a Array2D type
let Aarr = A.ToArray2D();;
val Aarr : float [,] = [[1.0; 7.0; 2.0]
[1.0; 3.0; 1.0]
[2.0; 9.0; 1.0]]
//Convert to vector (when the matrix has one row or column)
let Avec = A.ToVector()
let Arvec = A.ToRowVector();;
val Avec : Vector = vector [|1.0; 1.0; 2.0|]
val Arvec : RowVector = rowvec [|1.0; 7.0; 2.0|]

//operator [,]
A.[2,2];;
A.[2,2] <- 100.0;; A.[2,2];;
> val it : float = 100.0
> val it : unit = ()
> val it : float = 100.0

//sub-matrix
A.[1..2,1..2];;
A.[1..2,1..2] <- matrix [[2.;3.]; [8.;9.;]];; A.[1..2,1..2];;
> val it : Matrix = matrix [[3.0; 1.0]
[9.0; 1.0]]
> val it : unit = ()
> val it : Matrix = matrix [[2.0; 3.0]
[8.0; 9.0]]
//columns, rows
A.Column 2;;
A.Row 2;;
A.Columns (1,2);;
A.Rows (1,2);;
> val it : Vector = vector [|2.0; 3.0; 9.0|]
> val it : RowVector = rowvec [|2.0; 8.0; 9.0|]
> val it : Matrix = matrix [[7.0; 2.0]
[2.0; 3.0]
[8.0; 9.0]]
> val it : Matrix = matrix [[1.0; 2.0; 3.0]
[2.0; 8.0; 9.0]]
// sum and prod of all elements
let Asum = Matrix.sum A;;
let Aprod = Matrix.prod A;;
> val Asum : float = 35.0
> val Aprod : float = 12096.0

// create a matrix with 1
let C = Matrix.create 5 5 1.0;;
> val C : matrix = matrix [[1.0; 1.0; 1.0; 1.0; 1.0]
[1.0; 1.0; 1.0; 1.0; 1.0]
[1.0; 1.0; 1.0; 1.0; 1.0]
[1.0; 1.0; 1.0; 1.0; 1.0]
[1.0; 1.0; 1.0; 1.0; 1.0]]

// create a matrix with a function
let table = Matrix.init 3 3 (fun i j ->; (float i + 1.) * (float j + 1.));;
> val table : matrix = matrix [[1.0; 2.0; 3.0]
[2.0; 4.0; 6.0]
[3.0; 6.0; 9.0]]
// Identity
let I = Matrix.identity 2;;val I : matrix = matrix [[1.0; 0.0]
[0.0; 1.0]]
// trace
let Atrace = Matrix.trace A;;
val Atrace : float = 12.0

//Map a function to the elements
let Asqr = Matrix.map (fun x -> x*x) A;;
> val Asqr : matrix = matrix [[1.0; 49.0; 4.0]
[1.0; 4.0; 9.0]
[4.0; 64.0; 81.0]]
//Sparse matrix
let D = Matrix.initSparse 3 3 [ (0,0,1.0); (1,1,2.0); (2,2,3.0); ];;
> val D : matrix = matrix [[1.0; 0.0; 0.0]
[0.0; 2.0; 0.0]
[0.0; 0.0; 3.0]]
D.IsSparse;;
> val it : bool = true

//Matrix Generic
type intmatrix = Matrix
let intmatrix ll = Matrix.Generic.ofSeq ll
let C = intmatrix [ [1;2]; [3;4] ]
let D = intmatrix [ [1;1]; [1;1] ]
Matrix.Generic.inplaceAdd C D;;
C+D;;
C * 10;;type intmatrix = Matrix
val intmatrix : seq<#seq<'b>> -> Matrix<'b>
val C : Matrix = matrix [[2; 3]
[4; 5]]
val D : Matrix = matrix [[1; 1]
[1; 1]]
> val it : Matrix = matrix [[3; 4]
[5; 6]]
> val it : Matrix = matrix [[20; 30]
[40; 50]]
//Linear Algebra - float

//Solve linear equations
let A1 =
[ [1; 2; 3];
[1; 4; 9];
[1; 8; 27] ]
let A2 = Matrix.ofSeq(Seq.map (Seq.map float) A1)
let b = vector [2.; 3.; 4.]
let x = LinearAlgebra.Float.linear_solve A b;;
val A1 : int list list = [[1; 2; 3]; [1; 4; 9]; [1; 8; 27]]
val A2 : matrix = matrix [[1.0; 2.0; 3.0]
[1.0; 4.0; 9.0]
[1.0; 8.0; 27.0]]
val b : vector = vector [|2.0; 3.0; 4.0|]
val x : vector = vector [|0.5; 1.0; -0.1666666667|]

//Determinant
LinearAlgebra.Float.determinant A2;;
> val it : float = -12.0

//Inverse
LinearAlgebra.Float.inverse A2;;> val it : matrix = matrix [[3.0; -2.5; 0.5]
[-1.5; 2.0; -0.5]
[0.3333333333; -0.5; 0.1666666667]]
//svd
LinearAlgebra.Float.svd A2;;
> val it : matrix * matrix * matrix =
(matrix [[0.1166955507; 0.6917676817; 0.7126286712]
[0.3269819117; 0.6507677484; -0.6852621157]
[0.9377979409; -0.3129837252; 0.1502538182]],
matrix [[30.0404326; 0.0; 0.0]
[0.0; 1.878075836; 0.0]
[0.0; 0.0; 0.2126972812]],
matrix [[0.0459872007; 0.5481949583; 0.835085304]
[0.301051; 0.7894977011; -0.5348473383]
[0.9524985421; -0.2759993977; 0.1287278512]])

//Linear Algebra – rational

//Solve linear equations
let A3 =
[ [1; 2; 3];
[1; 4; 9];
[1; 8; 27] ]
let A4 = Matrix.Generic.ofSeq(Seq.map (Seq.map BigNum.FromInt) A3)
let b1 = Vector.Generic.ofSeq [2N; 3N; 4N]
let x1 = LinearAlgebra.Rational.linear_solve A4 b1;;
A4 * x1;;
> val A3 : int list list = [[1; 2; 3]; [1; 4; 9]; [1; 8; 27]]
val A4 : Matrix = matrix [[1N; 2N; 3N]
[1N; 4N; 9N]
[1N; 8N; 27N]]
val b1 : Vector = vector [|2N; 3N; 4N|]
val x1 : Vector = vector [|1/2N; 1N; -1/6N|]
> val it : Vector = vector [|2N; 3N; 4N|]


//Inverse - rational
LinearAlgebra.Rational.inverse A4;;
> val it : Matrix = matrix [[3N; -5/2N; 1/2N]
[-3/2N; 2N; -1/2N]
[1/3N; -1/2N; 1/6N]]

//Linear Algebra - complex

type cmatrix = Matrix
let cmatrix ll = Matrix.Generic.ofSeq ll

let A5 = cmatrix [ [ complex 1.0 1.0; complex 1.0 1.0];
[ complex 1.0 1.0; complex 1.0 1.0]];;
> type cmatrix = Matrix
val cmatrix : seq<#seq<'b>> -> Matrix<'b>
val A5 : Matrix = matrix [[1r+1i; 1r+1i]
[1r+1i; 1r+1i]]
//Determinant - complex
LinearAlgebra.Complex.determinant A5;;
> val it : Complex = 0r+0i {Conjugate = 0r+0i;
ImaginaryPart = 0.0;
Magnitude = 0.0;
Phase = 0.0;
RealPart = 0.0;
i = 0.0;
r = 0.0;}
f (A-B);;