3 people like it.

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

shows a simple implementation of a vector and matrix type together with a QR-decomposition using the Gram-Schmidt method. The algorithms themselfes are rather easy but I think the implementation of the types and the computations using recursive techniques might be interessting

Helpers

1: 
2: 
3: 
module Helpers =
    let curry f a b = f (a,b)
    let uncurry f (a,b) = f a b

Vector-type

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
 9: 
10: 
11: 
12: 
13: 
14: 
15: 
16: 
17: 
18: 
19: 
20: 
21: 
22: 
23: 
24: 
25: 
26: 
27: 
28: 
29: 
30: 
31: 
/// 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

Matrix - type

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
 9: 
10: 
11: 
12: 
13: 
14: 
15: 
16: 
17: 
18: 
19: 
20: 
21: 
22: 
23: 
24: 
25: 
26: 
27: 
28: 
29: 
30: 
31: 
32: 
33: 
34: 
35: 
36: 
37: 
38: 
39: 
40: 
41: 
42: 
43: 
44: 
45: 
46: 
47: 
48: 
49: 
50: 
51: 
52: 
53: 
54: 
55: 
56: 
57: 
58: 
59: 
60: 
61: 
62: 
63: 
64: 
65: 
/// 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 

Gram-Schmidt method

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
 9: 
10: 
11: 
12: 
13: 
14: 
15: 
16: 
17: 
18: 
19: 
20: 
21: 
22: 
23: 
24: 
25: 
26: 
/// 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

QR-decomposition

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
 9: 
10: 
11: 
12: 
13: 
/// 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
val curry : f:('a * 'b -> 'c) -> a:'a -> b:'b -> 'c

Full name: Script.Helpers.curry
val f : ('a * 'b -> 'c)
val a : 'a
val b : 'b
val uncurry : f:('a -> 'b -> 'c) -> a:'a * b:'b -> 'c

Full name: Script.Helpers.uncurry
val f : ('a -> 'b -> 'c)
Multiple items
union case Vector.Vector: float array -> Vector

--------------------
type Vector =
  | Vector of float array
  member Item : i:int -> float
  member Dim : int
  member Length : float
  member Normal : Vector
  member private Values : float array
  static member InnerProd : a:Vector * b:Vector -> float
  static member Len2 : a:Vector -> float
  static member Normalize : a:Vector -> Vector
  static member Zero : n:int -> Vector
  static member Len : (Vector -> float)
  static member ( + ) : a:Vector * b:Vector -> Vector
  static member ( .*. ) : a:Vector * b:Vector -> float
  static member ( * ) : s:float * a:Vector -> Vector
  static member ( - ) : a:Vector * b:Vector -> Vector

Full name: Script.Vector


 Represent a vector as an float array
Multiple items
val float : value:'T -> float (requires member op_Explicit)

Full name: Microsoft.FSharp.Core.Operators.float

--------------------
type float = System.Double

Full name: Microsoft.FSharp.Core.float

--------------------
type float<'Measure> = float

Full name: Microsoft.FSharp.Core.float<_>
type 'T array = 'T []

Full name: Microsoft.FSharp.Core.array<_>
val v : Vector
member private Vector.Values : float array

Full name: Script.Vector.Values


 shortcut to access the float array
val v' : float array
member Vector.Item : i:int -> float

Full name: Script.Vector.Item


 access a component of the vector
val i : int
property Vector.Values: float array


 shortcut to access the float array
member Vector.Dim : int

Full name: Script.Vector.Dim


 the vectors dimension = nr of components
module Array

from Microsoft.FSharp.Collections
val length : array:'T [] -> int

Full name: Microsoft.FSharp.Collections.Array.length
val a : Vector
val b : Vector
val zip : array1:'T1 [] -> array2:'T2 [] -> ('T1 * 'T2) []

Full name: Microsoft.FSharp.Collections.Array.zip
val map : mapping:('T -> 'U) -> array:'T [] -> 'U []

Full name: Microsoft.FSharp.Collections.Array.map
module Helpers

from Script
val s : float
static member Vector.InnerProd : a:Vector * b:Vector -> float

Full name: Script.Vector.InnerProd


 computes the inner product (also known as scalar product) of two vectors
val sumBy : projection:('T -> 'U) -> array:'T [] -> 'U (requires member ( + ) and member get_Zero)

Full name: Microsoft.FSharp.Collections.Array.sumBy
static member Vector.InnerProd : a:Vector * b:Vector -> float


 computes the inner product (also known as scalar product) of two vectors
static member Vector.Len2 : a:Vector -> float

Full name: Script.Vector.Len2


 the squared length of a vector
static member Vector.Len : (Vector -> float)

Full name: Script.Vector.Len


 the length of a vector
namespace System
type Math =
  static val PI : float
  static val E : float
  static member Abs : value:sbyte -> sbyte + 6 overloads
  static member Acos : d:float -> float
  static member Asin : d:float -> float
  static member Atan : d:float -> float
  static member Atan2 : y:float * x:float -> float
  static member BigMul : a:int * b:int -> int64
  static member Ceiling : d:decimal -> decimal + 1 overload
  static member Cos : d:float -> float
  ...

Full name: System.Math
System.Math.Sqrt(d: float) : float
static member Vector.Len2 : a:Vector -> float


 the squared length of a vector
member Vector.Length : float

Full name: Script.Vector.Length
property Vector.Len: Vector -> float


 the length of a vector
static member Vector.Normalize : a:Vector -> Vector

Full name: Script.Vector.Normalize


 normalize a vector (the result will have the same direction but length 1)
property Vector.Length: float
member Vector.Normal : Vector

Full name: Script.Vector.Normal
static member Vector.Normalize : a:Vector -> Vector


 normalize a vector (the result will have the same direction but length 1)
static member Vector.Zero : n:int -> Vector

Full name: Script.Vector.Zero


 create a zero-vector
val n : int
val zeroCreate : count:int -> 'T []

Full name: Microsoft.FSharp.Collections.Array.zeroCreate
Multiple items
union case Matrix.Matrix: float [,] -> Matrix

--------------------
type Matrix =
  | Matrix of float [,]
  member Column : i:int -> Vector
  member Columns : Vector []
  member NrColumns : int
  member NrRows : int
  member private Values : float [,]
  static member ColVectors : m:Matrix -> Vector []
  static member Create : vs:Vector array -> Matrix
  static member Create : rows:int * cols:int * f:(int -> int -> float) -> Matrix
  static member Transpose : m:Matrix -> Matrix
  static member Zero : rows:int * cols:int -> 'a [,]
  static member ( + ) : a:Matrix * b:Matrix -> Matrix
  static member ( * ) : a:Matrix * b:Matrix -> Matrix
  static member ( * ) : s:float * b:Matrix -> Matrix
  static member ( - ) : a:Matrix * b:Matrix -> Matrix

Full name: Script.Matrix


 represent a matrix as a 2-dimensional float array
val m : Matrix
member private Matrix.Values : float [,]

Full name: Script.Matrix.Values


 shortcut to the array
val m' : float [,]
member Matrix.NrRows : int

Full name: Script.Matrix.NrRows


 the number of rows
module Array2D

from Microsoft.FSharp.Collections
val length1 : array:'T [,] -> int

Full name: Microsoft.FSharp.Collections.Array2D.length1
property Matrix.Values: float [,]


 shortcut to the array
member Matrix.NrColumns : int

Full name: Script.Matrix.NrColumns


 the number of columns
val length2 : array:'T [,] -> int

Full name: Microsoft.FSharp.Collections.Array2D.length2
member Matrix.Column : i:int -> Vector

Full name: Script.Matrix.Column


 access a column-vector via its zero-based index
val z : int
member Matrix.Columns : Vector []

Full name: Script.Matrix.Columns


 get all column-vectors as an array
member Matrix.Column : i:int -> Vector


 access a column-vector via its zero-based index
static member Matrix.ColVectors : m:Matrix -> Vector []

Full name: Script.Matrix.ColVectors
property Matrix.Columns: Vector []


 get all column-vectors as an array
val a : Matrix
val b : Matrix
val init : length1:int -> length2:int -> initializer:(int -> int -> 'T) -> 'T [,]

Full name: Microsoft.FSharp.Collections.Array2D.init
property Matrix.NrRows: int


 the number of rows
property Matrix.NrColumns: int


 the number of columns
val r : int
val c : int
Multiple items
module List

from Microsoft.FSharp.Collections

--------------------
type List<'T> =
  | ( [] )
  | ( :: ) of Head: 'T * Tail: 'T list
  interface IEnumerable
  interface IEnumerable<'T>
  member Head : 'T
  member IsEmpty : bool
  member Item : index:int -> 'T with get
  member Length : int
  member Tail : 'T list
  static member Cons : head:'T * tail:'T list -> 'T list
  static member Empty : 'T list

Full name: Microsoft.FSharp.Collections.List<_>
val map : mapping:('T -> 'U) -> list:'T list -> 'U list

Full name: Microsoft.FSharp.Collections.List.map
val sum : list:'T list -> 'T (requires member ( + ) and member get_Zero)

Full name: Microsoft.FSharp.Collections.List.sum
static member Matrix.Transpose : m:Matrix -> Matrix

Full name: Script.Matrix.Transpose


 computes the transpose of a matrix
val x : int
val y : int
static member Matrix.Create : vs:Vector array -> Matrix

Full name: Script.Matrix.Create


 creates a matrix from a array of columnvectors
val vs : Vector array
val l : int
property System.Array.Length: int
static member Matrix.Create : rows:int * cols:int * f:(int -> int -> float) -> Matrix

Full name: Script.Matrix.Create


 creates a matrix with an generator-function (the row r and col c of the matrix will have a value of f(r,c))
val rows : int
val cols : int
val f : (int -> int -> float)
static member Matrix.Zero : rows:int * cols:int -> 'a [,]

Full name: Script.Matrix.Zero


 creates a zero-Matrix
val zeroCreate : length1:int -> length2:int -> 'T [,]

Full name: Microsoft.FSharp.Collections.Array2D.zeroCreate
val private project : baseV:Vector -> v:Vector -> Vector

Full name: Script.GramSchmidt.project


 projects a vector onto another vector
val baseV : Vector
val Orthogonalize : vs:Vector list -> Vector list

Full name: Script.GramSchmidt.Orthogonalize


 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
val vs : Vector list
type 'T list = List<'T>

Full name: Microsoft.FSharp.Collections.list<_>
val calcNextVec : (Vector -> Vector list -> Vector)
val cur : Vector
val found : Vector list
val vs' : Vector list
val calcVecs : (Vector list -> Vector list -> Vector list)
val toDo : Vector list
val rest : Vector list
val rev : list:'T list -> 'T list

Full name: Microsoft.FSharp.Collections.List.rev
val OrthogonalizeNormal : (Vector list -> Vector list)

Full name: Script.GramSchmidt.OrthogonalizeNormal


 computes a set of orthonormal vectors that span the same subspace as the vectors in vs
module QRdecomposition

from Script


 decomposes a square matrix into a product of an orthongal and an upper triangular matrix
val UsingGramSchmidt : m:Matrix -> Matrix * Matrix

Full name: Script.QRdecomposition.UsingGramSchmidt


 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
val colVs : Vector []
val normVs : Vector []
val ofList : list:'T list -> 'T []

Full name: Microsoft.FSharp.Collections.Array.ofList
module GramSchmidt

from Script


 computes a set of orthogonal or orthonormal vectors that spans the same
 subspace as the given vector-set
val ofArray : array:'T [] -> 'T list

Full name: Microsoft.FSharp.Collections.List.ofArray
val Q : Matrix
static member Matrix.Create : vs:Vector array -> Matrix


 creates a matrix from a array of columnvectors

static member Matrix.Create : rows:int * cols:int * f:(int -> int -> float) -> 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))
val R : Matrix
Raw view Test code New version

More information

Link:http://fssnip.net/3p
Posted:13 years ago
Author:Carsten König
Tags: qr-decomposition , gram-schmidt , linear algebra , matrix , vector