Jeremy Bellows - Gram Schmidt in FSharp

Gram Schmidt in FSharp

-


Gram Schmidt

The Gram-Schmidt process is a function that accepts a set of vectors and returns the corresponding orthonormal vectors. Acceptable sets of vectors are finite, linearly independent, and identical in vector sizes.

Gram Schmidt in F

Below is the implementation of the Gram-Schmidt process in F#.
Feedback is encourage and appreciated!
Tweet @JeremyBellows
Leave a comment or star the Github Gist of this code

type Vector = {Vector : float list} with
  static member Create vector =
    { Vector = vector }

  member this.InnerProduct otherVector =
    List.zip this.Vector otherVector.Vector
    |> List.map (fun (vectorOneNumber, vectorTwoNumber) -> vectorOneNumber * vectorTwoNumber)
    |> List.sum

  member this.Length =
    let vectorLengthSummedWithSquaredValues =
      this.Vector
      |> List.map (fun scalar -> pown scalar 2)
      |> List.sum
    sqrt vectorLengthSummedWithSquaredValues

  member this.MultiplyByScalar scalar =
    this.Vector
    |> List.map (fun vectorNumber -> vectorNumber * scalar)
    |> Vector.Create

  member this.DivideByScalar scalar =
    this.Vector
    |> List.map (fun vectorScalar -> vectorScalar / scalar)
    |> Vector.Create

  member this.CalculateProjection vector =
    let innerProductScalar =
      let innerProductOfProjectionVectorAgainstVector = this.InnerProduct vector
      let innerProductOfVector = vector.InnerProduct vector
      innerProductOfProjectionVectorAgainstVector / innerProductOfVector
    vector.MultiplyByScalar innerProductScalar

  member this.SubtractVector vector =
    let subtractNumbers numberOne numberTwo =
      numberOne - numberTwo
    List.map2 subtractNumbers this.Vector vector.Vector
    |> Vector.Create

type ListOfVectors = Vector list

let calculateOrthogonalVectors (vectors : ListOfVectors) : ListOfVectors =
  let rec calculateOrthogonalVectors (vectors : ListOfVectors) (orthogonalVectors : ListOfVectors) =
    let findOrthagonalVector (vector : Vector) (previousOrthogonalVectors : ListOfVectors) : Vector =
      let projections =
        previousOrthogonalVectors |> List.map vector.CalculateProjection
      let calculateOrthogonalVector (accumulatorVector : Vector) (projection : Vector) =
        accumulatorVector.SubtractVector projection
      projections |> List.fold calculateOrthogonalVector vector
    match vectors with
      | [] ->
        orthogonalVectors
      | vector :: remainingVectors ->
        let updatedOrthogonalVectors : ListOfVectors =
          if orthogonalVectors |> List.isEmpty then
            //Vector is orthogonal to itself
            [vector]
          else
            let orthagonalVector = findOrthagonalVector vector orthogonalVectors
            List.append orthogonalVectors [orthagonalVector]
        calculateOrthogonalVectors remainingVectors updatedOrthogonalVectors
  calculateOrthogonalVectors vectors []

let calculateOrthonormalVectors (orthogonalVectors : ListOfVectors) : ListOfVectors =
  let rec calculateOrthonormalVectors (orthogonalVectors : ListOfVectors) orthonormalVectors =
    match orthogonalVectors with
     | [] ->
      orthonormalVectors
     | orthogonalVector :: remainingOrthogonalVectors ->
      let updatedOrthonormalVectors =
        let orthonormalVector =
          orthogonalVector.DivideByScalar orthogonalVector.Length
        List.append orthonormalVectors [orthonormalVector]
      calculateOrthonormalVectors remainingOrthogonalVectors updatedOrthonormalVectors
  calculateOrthonormalVectors orthogonalVectors []

let (==) floatNumber (integerNumber : int) =
  System.BitConverter.DoubleToInt64Bits(floatNumber) = (int64 integerNumber)

let isOrthonormal (vectors : ListOfVectors) =
  let checkMagnitude (vector : Vector) =
    vector.Length == 1
  let vectorsHaveMagnitudeOfOne = vectors |> List.forall checkMagnitude
  let rec checkDotProducts (remainingVectors : ListOfVectors) =
    match remainingVectors with
    | [] -> true
    | vector :: remainingVectors ->
      let checkDotProduct comparisonVector =
        (vector.InnerProduct comparisonVector) == 0
      match remainingVectors |> List.forall checkDotProduct with
      | true -> checkDotProducts remainingVectors
      | false -> false
  let vectorsHaveZeroDotProducts = vectors |> checkDotProducts
  vectorsHaveMagnitudeOfOne && vectorsHaveZeroDotProducts

//Code to test sanity
let convertToVectorList floatLists =
  floatLists
  |> List.map Vector.Create

let brokenVector = [[1.999; 0.0; 0.0]; [0.0; 1.999; 0.0]] |> convertToVectorList
isOrthonormal brokenVector

let testVectorsOne =
  [[1.0; 1.0; 1.0; 1.0]; [1.0; 0.0; 0.0; 1.0]; [0.0; 2.0; 1.0; -1.0]]
  |> convertToVectorList
let testOneOrthogonalVectors = calculateOrthogonalVectors testVectorsOne
let testOneOrthonormalVectors = calculateOrthonormalVectors testOneOrthogonalVectors
isOrthonormal testOneOrthonormalVectors

let testVectorsTwo =
  [[1.0; 0.0; 0.0]; [2.0; 0.0; 3.0]; [4.0;5.0;6.0]]
  |> convertToVectorList
let testTwoOrthogonalVectors = calculateOrthogonalVectors testVectorsTwo
let testTwoOrthonormalVectors = calculateOrthonormalVectors testTwoOrthogonalVectors
isOrthonormal testTwoOrthonormalVectors

let testVectorsThree =
  [[2.0; 2.0; 1.0]; [-2.0; 1.0; 2.0]; [18.0; 0.0; 0.0]]
  |> convertToVectorList
let testThreeOrthogonalVectors = calculateOrthogonalVectors testVectorsThree
let testThreeOrthonormalVectors = calculateOrthonormalVectors testThreeOrthogonalVectors
isOrthonormal testThreeOrthonormalVectors

let testVectorsFour =
  [[2.0; 2.0; 1.0]; [1.0; 1.0; 5.0]]
  |> convertToVectorList
let testFourOrthogonalVectors = calculateOrthogonalVectors testVectorsFour
let testFourOrthonormalVectors = calculateOrthonormalVectors testFourOrthogonalVectors
isOrthonormal testFourOrthonormalVectors