How to resolve the algorithm Cholesky decomposition step by step in the F# programming language

Published on 12 May 2024 09:40 PM
#F_Sharp

How to resolve the algorithm Cholesky decomposition step by step in the F# programming language

Table of Contents

Problem Statement

Every symmetric, positive definite matrix A can be decomposed into a product of a unique lower triangular matrix L and its transpose:

L

{\displaystyle L}

is called the Cholesky factor of

A

{\displaystyle A}

, and can be interpreted as a generalized square root of

A

{\displaystyle A}

, as described in Cholesky decomposition. In a 3x3 example, we have to solve the following system of equations: We can see that for the diagonal elements (

l

k k

{\displaystyle l_{kk}}

) of

L

{\displaystyle L}

there is a calculation pattern: or in general: For the elements below the diagonal (

l

i k

{\displaystyle l_{ik}}

, where

i

k

{\displaystyle i>k}

) there is also a calculation pattern: which can also be expressed in a general formula: Task description The task is to implement a routine which will return a lower Cholesky factor

L

{\displaystyle L}

for every given symmetric, positive definite nxn matrix

A

{\displaystyle A}

. You should then test it on the following two examples and include your output. Example 1: Example 2:

Let's start with the solution:

Step by Step solution about How to resolve the algorithm Cholesky decomposition step by step in the F# programming language

Source code in the fsharp programming language

open Microsoft.FSharp.Collections

let cholesky a =
    let calc (a: float[,]) (l: float[,]) i j =
        let c1 j =
            let sum = List.sumBy (fun k -> l.[j, k] ** 2.0) [0..j - 1]
            sqrt (a.[j, j] - sum)
        let c2 i j = 
            let sum = List.sumBy (fun k -> l.[i, k] * l.[j, k]) [0..j - 1]
            (1.0 / l.[j, j]) * (a.[i, j] - sum)
        if j > i then 0.0 else
            if i = j
            then c1 j
            else c2 i j
    let l = Array2D.zeroCreate (Array2D.length1 a) (Array2D.length2 a)
    Array2D.iteri (fun i j _ -> l.[i, j] <- calc a l i j) l
    l

let printMat a =
    let arrow = (Array2D.length2 a |> float) / 2.0 |> int
    let c = cholesky a
    for row in 0..(Array2D.length1 a) - 1 do
        for col in 0..(Array2D.length2 a) - 1 do
            printf "%.5f,\t" a.[row, col]
        printf (if arrow = row then "--> \t" else "\t\t")
        for col in 0..(Array2D.length2 c) - 1 do
            printf "%.5f,\t" c.[row, col]
        printfn ""

let ex1 = array2D [
    [25.0; 15.0; -5.0];
    [15.0; 18.0; 0.0];
    [-5.0; 0.0; 11.0]]

let ex2 = array2D [
    [18.0; 22.0; 54.0; 42.0];
    [22.0; 70.0; 86.0; 62.0];
    [54.0; 86.0; 174.0; 134.0];
    [42.0; 62.0; 134.0; 106.0]]

printfn "ex1:"
printMat ex1

printfn "ex2:"
printMat ex2

You may also check:How to resolve the algorithm Matrix multiplication step by step in the Liberty BASIC programming language
You may also check:How to resolve the algorithm Loops/Wrong ranges step by step in the Racket programming language
You may also check:How to resolve the algorithm String case step by step in the ALGOL W programming language
You may also check:How to resolve the algorithm Window creation step by step in the Nim programming language
You may also check:How to resolve the algorithm Rename a file step by step in the Plain English programming language