Cholesky decomposition
You are encouraged to solve this task according to the task description, using any language you may know.
Method description
Every symmetric, positive definite matrix A can be decomposed into a product of a unique lower triangular matrix L and its transpose:
is called the Cholesky factor of , and can be interpreted as a generalized square root of , 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 () of there is a calculation pattern:
or in general:
For the elements below the diagonal (, where ) 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 for every given symmetric, positive definite nxn matrix . You should then test it on the following two examples and include your output.
Example 1:
25 15 -5 5 0 0 15 18 0 --> 3 3 0 -5 0 11 -1 1 3
Example 2:
18 22 54 42 4.24264 0.00000 0.00000 0.00000 22 70 86 62 --> 5.18545 6.56591 0.00000 0.00000 54 86 174 134 12.72792 3.04604 1.64974 0.00000 42 62 134 106 9.89949 1.62455 1.84971 1.39262
Common Lisp
<lang lisp>;; 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))</lang>
<lang lisp>
- 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))
</lang>
<lang lisp>
- 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))
</lang>