How to resolve the algorithm Cholesky decomposition step by step in the Common Lisp programming language
How to resolve the algorithm Cholesky decomposition step by step in the Common Lisp 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 Common Lisp programming language
Source code in the common programming language
;; Calculates the Cholesky decomposition matrix L
;; for a positive-definite, symmetric nxn matrix A.
(defun chol (A)
(let* ((n (car (array-dimensions A)))
(L (make-array `(,n ,n) :initial-element 0)))
(do ((k 0 (incf k))) ((> k (- n 1)) nil)
;; First, calculate diagonal elements L_kk.
(setf (aref L k k)
(sqrt (- (aref A k k)
(do* ((j 0 (incf j))
(sum (expt (aref L k j) 2)
(incf sum (expt (aref L k j) 2))))
((> j (- k 1)) sum)))))
;; Then, all elements below a diagonal element, L_ik, i=k+1..n.
(do ((i (+ k 1) (incf i)))
((> i (- n 1)) nil)
(setf (aref L i k)
(/ (- (aref A i k)
(do* ((j 0 (incf j))
(sum (* (aref L i j) (aref L k j))
(incf sum (* (aref L i j) (aref L k j)))))
((> j (- k 1)) sum)))
(aref L k k)))))
;; Return the calculated matrix L.
L))
;; Example 1:
(setf A (make-array '(3 3) :initial-contents '((25 15 -5) (15 18 0) (-5 0 11))))
(chol A)
#2A((5.0 0 0)
(3.0 3.0 0)
(-1.0 1.0 3.0))
;; Example 2:
(setf B (make-array '(4 4) :initial-contents '((18 22 54 42) (22 70 86 62) (54 86 174 134) (42 62 134 106))))
(chol B)
#2A((4.2426405 0 0 0)
(5.18545 6.565905 0 0)
(12.727922 3.0460374 1.6497375 0)
(9.899495 1.6245536 1.849715 1.3926151))
;; case of matrix stored as a list of lists (inner lists are rows of matrix)
;; as above, returns the Cholesky decomposition matrix of a square positive-definite, symmetric matrix
(defun cholesky (m)
(let ((l (list (list (sqrt (caar m))))) x (j 0) i)
(dolist (cm (cdr m) (mapcar #'(lambda (x) (nconc x (make-list (- (length m) (length x)) :initial-element 0))) l))
(setq x (list (/ (car cm) (caar l))) i 0)
(dolist (cl (cdr l))
(setf (cdr (last x)) (list (/ (- (elt cm (incf i)) (*v x cl)) (car (last cl))))))
(setf (cdr (last l)) (list (nconc x (list (sqrt (- (elt cm (incf j)) (*v x x))))))))))
;; where *v is the scalar product defined as
(defun *v (v1 v2) (reduce #'+ (mapcar #'* v1 v2)))
;; example 1
CL-USER> (setf a '((25 15 -5) (15 18 0) (-5 0 11)))
((25 15 -5) (15 18 0) (-5 0 11))
CL-USER> (cholesky a)
((5 0 0) (3 3 0) (-1 1 3))
CL-USER> (format t "~{~{~5d~}~%~}" (cholesky a))
5 0 0
3 3 0
-1 1 3
NIL
;; example 2
CL-USER> (setf a '((18 22 54 42) (22 70 86 62) (54 86 174 134) (42 62 134 106)))
((18 22 54 42) (22 70 86 62) (54 86 174 134) (42 62 134 106))
CL-USER> (cholesky a)
((4.2426405 0 0 0) (5.18545 6.565905 0 0) (12.727922 3.0460374 1.6497375 0) (9.899495 1.6245536 1.849715 1.3926151))
CL-USER> (format t "~{~{~10,5f~}~%~}" (cholesky a))
4.24264 0.00000 0.00000 0.00000
5.18545 6.56591 0.00000 0.00000
12.72792 3.04604 1.64974 0.00000
9.89950 1.62455 1.84971 1.39262
NIL
You may also check:How to resolve the algorithm Sort an outline at every level step by step in the Wren programming language
You may also check:How to resolve the algorithm State name puzzle step by step in the C programming language
You may also check:How to resolve the algorithm Run-length encoding step by step in the PureBasic programming language
You may also check:How to resolve the algorithm Hash from two arrays step by step in the EchoLisp programming language
You may also check:How to resolve the algorithm Catalan numbers step by step in the Logtalk programming language