QR-decomoposition of a square-matrix using the Gram-Schmidt method

module Helpers =
    let curry f a b = f (a,b)
    let uncurry f (a,b) = f a b

/// Represent a vector as an float array
type Vector =  | Vector of float array
    with 
    /// shortcut to access the float array
    member private v.Values = match v with Vector v' -> v'
    /// access a component of the vector
    member v.Item(i) = v.Values.[i]
    /// the vectors dimension = nr of components
    member v.Dim = Array.length v.Values

    // your normal vector-operators
    static member (+) (a : Vector, b : Vector) = Array.zip a.Values b.Values |> Array.map (Helpers.uncurry (+)) |> Vector
    static member (-) (a : Vector, b : Vector) = Array.zip a.Values b.Values |> Array.map (Helpers.uncurry (-)) |> Vector
    static member (*) (s : float, a : Vector) = a.Values |> Array.map ((*) s) |> Vector

    /// computes the inner product (also known as scalar product) of two vectors
    static member InnerProd (a : Vector, b : Vector) = Array.zip a.Values b.Values |> Array.sumBy (Helpers.uncurry (*))
    static member (.*.) (a,b) = Vector.InnerProd(a,b)

    /// the squared length of a vector
    static member Len2 (a : Vector) = a .*. a
    /// the length of a vector
    static member Len = System.Math.Sqrt << Vector.Len2
    member v.Length = Vector.Len v

    /// normalize a vector (the result will have the same direction but length 1)
    static member Normalize (a : Vector) = (1.0 / a.Length) * a
    member v.Normal = Vector.Normalize v

    /// create a zero-vector
    static member Zero(n) = Array.zeroCreate n |> Vector

/// represent a matrix as a 2-dimensional float array
type Matrix = Matrix of float [,]
    with
    /// shortcut to the array
    member private m.Values = match m  with Matrix m' -> m'
    
    /// the number of rows
    member m.NrRows = Array2D.length1 m.Values
    /// the number of columns
    member m.NrColumns = Array2D.length2 m.Values

    /// access a column-vector via its zero-based index
    member m.Column(i) = [| 0..Array2D.length1 m.Values-1 |] |> Array.map (fun z -> m.Values.[z,i]) |> Vector

    /// get all column-vectors as an array
    member m.Columns = [|0..Array2D.length2 m.Values-1 |] |> Array.map m.Column
    static member ColVectors (m : Matrix) = m.Columns

    /// Matrix-multiplication
    static member (*) (a : Matrix, b : Matrix) =
        Array2D.init 
            a.NrRows b.NrColumns
            (fun r c -> [0..a.NrColumns-1] 
                        |> List.map (fun i -> a.Values.[r,i] * b.Values.[i,c])
                        |> List.sum)
        |> Matrix
        
    /// Matrix-addition
    static member (+) (a : Matrix, b : Matrix) =
        Array2D.init 
            a.NrRows b.NrColumns
            (fun r c -> a.Values.[r,c] + b.Values.[r,c])
        |> Matrix
        
    /// Matrix-substraction
    static member (-) (a : Matrix, b : Matrix) =
        Array2D.init 
            a.NrRows b.NrColumns
            (fun r c -> a.Values.[r,c] + b.Values.[r,c])
        |> Matrix
        
    /// scalar - multiplication
    static member (*) (s : float, b : Matrix) =
        Array2D.init 
            b.NrRows b.NrColumns
            (fun r c -> s * b.Values.[r,c])
        |> Matrix
        

    /// computes the transpose of a matrix
    static member Transpose (m : Matrix) = Array2D.init m.NrRows m.NrColumns (fun x y -> m.Values.[y,x]) |> Matrix

    /// creates a matrix from a array of columnvectors
    static member Create (vs : Vector array) =
        let l = vs.[0].Dim
        Array2D.init l vs.Length (fun r c -> vs.[c].[r])
        |> Matrix

    /// creates a matrix with an generator-function (the row r and col c of the matrix will have a value of f(r,c))
    static member Create (rows, cols, f) =
        Array2D.init rows cols f
        |> Matrix

    /// creates a zero-Matrix
    static member Zero (rows, cols) = Array2D.zeroCreate rows cols 

/// computes a set of orthogonal or orthonormal vectors that spans the same
/// subspace as the given vector-set
module GramSchmidt =

    /// projects a vector onto another vector
    let private project (baseV : Vector) (v : Vector) : Vector 
        = (v .*. baseV)/(baseV .*. baseV) * baseV

    /// computes a set of orthogonal vectors that span the same subspace as the vectors in vs
    /// see the Wikipedia-article at http://en.wikipedia.org/wiki/Gram_schmidt for a good
    /// explanation of this method
    let Orthogonalize (vs : Vector list) =
        // computes one vector based on the allready found vectors
        let rec calcNextVec cur found =
            match found with
            | [] -> cur
            | v::vs' -> calcNextVec (cur - project v cur) vs'
        // computes the vectors
        let rec calcVecs toDo found =
            match toDo with
            | [] -> found
            | cur::rest -> calcVecs rest ((calcNextVec cur found)::found)
        calcVecs vs [] |> List.rev

    /// computes a set of orthonormal vectors that span the same subspace as the vectors in vs
    let OrthogonalizeNormal = Orthogonalize >> List.map Vector.Normalize

/// decomposes a square matrix into a product of an orthongal and an upper triangular matrix
module QRdecomposition =
    
    /// computes the QR-decomposition using the Gram-Schmidt method
    /// see the Wikipedia-article at http://en.wikipedia.org/wiki/QR_decomposition for a good
    /// explanation of this method
    let UsingGramSchmidt (m : Matrix) =
        let colVs = m.Columns
        let normVs = Array.ofList <| GramSchmidt.OrthogonalizeNormal (List.ofArray colVs)
        let Q = Matrix.Create normVs
        let f r c = if r <= c then normVs.[r] .*. colVs.[c] else 0.0
        let R = Matrix.Create (m.NrRows, m.NrColumns, f)
        Q,R