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

Published on 7 June 2024 03:52 AM
#Haskell

How to resolve the algorithm Cholesky decomposition step by step in the Haskell 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 Haskell programming language

The provided code in Haskell implements the Cholesky decomposition of a real, symmetric, and positive-definite matrix. The Cholesky decomposition is a factorization of a matrix into a lower triangular matrix and its transpose. This decomposition is useful for solving systems of linear equations and other matrix computations.

Data Structures:

  • Arr: Represents a two-dimensional unboxed array, where the indices are tuples of integers (Int, Int), and the elements are of type Double. This array represents the matrix to be decomposed.
  • Idx: Type alias for (Int, Int), representing indices in the matrix.

Main Function:

  • cholesky: Takes an Arr and computes its Cholesky decomposition, returning a lower triangular matrix as an Arr.

Helper Functions:

  • get: Retrieves an element from the lower triangular matrix.
  • maxBnd: Finds the maximum bounds of the input matrix.
  • update: Updates an element in the lower triangular matrix.

Example Matrices:

  • ex1 and ex2: Example matrices used for testing.

Main Function:

  • main:
    • For ex1 and ex2, it prints the lower triangular matrix resulting from the Cholesky decomposition.
    • For matrices a and b, it prints the original matrix, its Cholesky decomposition, and the transpose of the Cholesky decomposition.

Additional Details:

  • The Data.Array libraries provide various array data structures and operations.
  • unsafeFreeze is used to freeze an STUArray into a UArray for performance reasons.
  • splitAt is used to split a list into two parts at a given index for matrix formatting.
  • showMatrice is a helper function to format the matrix for printing.
  • The Numeric.LinearAlgebra library is also used for matrix operations.

Source code in the haskell programming language

module Cholesky (Arr, cholesky) where

import Data.Array.IArray
import Data.Array.MArray
import Data.Array.Unboxed
import Data.Array.ST

type Idx = (Int,Int)
type Arr = UArray Idx Double

-- Return the (i,j) element of the lower triangular matrix.  (We assume the
-- lower array bound is (0,0).)
get :: Arr -> Arr -> Idx -> Double
get a l (i,j) | i == j = sqrt $ a!(j,j) - dot
              | i  > j = (a!(i,j) - dot) / l!(j,j)
              | otherwise = 0
  where dot = sum [l!(i,k) * l!(j,k) | k <- [0..j-1]]

-- Return the lower triangular matrix of a Cholesky decomposition.  We assume
-- the input is a real, symmetric, positive-definite matrix, with lower array
-- bounds of (0,0).
cholesky :: Arr -> Arr
cholesky a = let n = maxBnd a
             in runSTUArray $ do
               l <- thaw a
               mapM_ (update a l) [(i,j) | i <- [0..n], j <- [0..n]]
               return l
  where maxBnd = fst . snd . bounds
        update a l i = unsafeFreeze l >>= \l' -> writeArray l i (get a l' i)


import Data.Array.IArray
import Data.List
import Cholesky

fm _ [] = "" 
fm _ [x] = fst x
fm width ((a,b):xs) = a ++ (take (width - b) $ cycle " ") ++ (fm width xs)

fmt width row (xs,[]) = fm width xs
fmt width row (xs,ys) = fm width xs  ++ "\n" ++ fmt width row (splitAt row ys)

showMatrice row xs   = ys where
  vs = map (\s -> let sh = show s in (sh,length sh)) xs
  width = (maximum $ snd $ unzip vs) + 1
  ys = fmt width row (splitAt row vs)

ex1, ex2 :: Arr
ex1 = listArray ((0,0),(2,2)) [25, 15, -5, 
                               15, 18,  0, 
                               -5,  0, 11]

ex2 = listArray ((0,0),(3,3)) [18, 22,  54,  42, 
                               22, 70,  86,  62, 
                               54, 86, 174, 134, 
                               42, 62, 134, 106]

main :: IO ()
main = do
  putStrLn $ showMatrice 3 $ elems $ cholesky ex1
  putStrLn $ showMatrice 4 $ elems $ cholesky ex2


import Numeric.LinearAlgebra

a,b :: Matrix R
a = (3><3) 
  [25, 15, -5
  ,15, 18, 0
  ,-5,  0, 11]

b = (4><4)
  [ 18, 22, 54, 42
  , 22, 70, 86, 62
  , 54, 86,174,134
  , 42, 62,134,106]

main = do
  let sa = sym a
      sb = sym b
  print sa
  print $ chol sa
  print sb
  print $ chol sb
  print $ tr $ chol sb

You may also check:How to resolve the algorithm History variables step by step in the Swift programming language
You may also check:How to resolve the algorithm HTTP step by step in the ActionScript programming language
You may also check:How to resolve the algorithm Idoneal numbers step by step in the XPL0 programming language
You may also check:How to resolve the algorithm 24 game step by step in the Logo programming language
You may also check:How to resolve the algorithm Formatted numeric output step by step in the ALGOL 68 programming language