### Author Topic: Challenge: Cholesky decomposition  (Read 14074 times)

0 Members and 1 Guest are viewing this topic.

#### ElpanovEvgeniy

• Water Moccasin
• Posts: 1569
• Moscow (Russia)
##### Challenge: Cholesky decomposition
« on: March 10, 2010, 09:59:32 AM »
New challenge!

You need - Cholesky decomposition of matrices.
Cholesky algorithm:
http://en.wikipedia.org/wiki/Cholesky_decomposition

Example of the original matrix:
Code: [Select]
`'((4.0 2.0 6.0 8.0 2.0)  (2.0 10.0 9.0 13.0 7.0)  (6.0 9.0 14.0 20.0 12.0)  (8.0 13.0 20.0 54.0 35.0)  (2.0 7.0 12.0 35.0 43.0) )`
result of the program:
Code: [Select]
`'((2.0 1.0 3.0 4.0 1.0)  (0.0 3.0 2.0 3.0 2.0)  (0.0 0.0 1.0 2.0 5.0)  (0.0 0.0 0.0 5.0 3.0)  (0.0 0.0 0.0 0.0 2.0) )`

#### Lee Mac

• Seagull
• Posts: 12943
• London, England
##### Re: Challenge: Cholesky decomposition
« Reply #1 on: March 10, 2010, 10:10:10 AM »
Great challenge Evgeniy

#### Swift

• Swamp Rat
• Posts: 596
##### Re: Challenge: Cholesky decomposition
« Reply #2 on: March 10, 2010, 10:17:34 AM »

#### pkohut

• Guest
##### Re: Challenge: Cholesky decomposition
« Reply #3 on: March 10, 2010, 10:23:35 AM »
Oh my, nice.

Can someone translate this for me? It's from the wiki link.
Quote
The algorithms described below all involve about n3/3 FLOPs, where n is the size of the matrix A
I don't have anything to reference it against.

#### Swift

• Swamp Rat
• Posts: 596
##### Re: Challenge: Cholesky decomposition
« Reply #4 on: March 10, 2010, 10:30:42 AM »
Oh my, nice.

Can someone translate this for me? It's from the wiki link.
Quote
The algorithms described below all involve about n3/3 FLOPs, where n is the size of the matrix A
I don't have anything to reference it against.

I think that is an estimate of the number of Floating Point Operations

#### Lee Mac

• Seagull
• Posts: 12943
• London, England
##### Re: Challenge: Cholesky decomposition
« Reply #5 on: March 10, 2010, 11:51:14 AM »
My first attempt:

Code: [Select]
`(defun LM_Cholesky (mat / v+v sxv put i j)  ;; Lee Mac  ~  10.03.10  (defun v-v (v1 v2)    (mapcar (function -) v1 v2))  (defun sxv (s v1)    (mapcar (function (lambda (x) (* s x))) v1))  (defun put (x n lst / k)    (setq k -1)    (mapcar (function (lambda (y) (if (= (setq k (1+ k)) n) x y))) lst))      (setq i -1)     (mapcar    (function      (lambda (row)        (setq i (1+ i) j -1)        (repeat i                    (setq mat (put (setq row (v-v row (sxv (/ (nth (setq j (1+ j)) row)                                                    (float (nth j (nth j mat))))                                                 (nth j mat))))                         i mat)          )                    )      )    )    mat) mat)          `
« Last Edit: March 10, 2010, 12:50:21 PM by Lee Mac »

#### ElpanovEvgeniy

• Water Moccasin
• Posts: 1569
• Moscow (Russia)
##### Re: Challenge: Cholesky decomposition
« Reply #6 on: March 10, 2010, 03:07:09 PM »
Hello Lee Mac
You turned interesting program, but there is a small mistake ...
For the test case:
Code: [Select]
`(LM_Cholesky '((4. 2. 6. 8. 2.) (2. 10. 9. 13. 7.) (6. 9. 14. 20. 12.) (8. 13. 20. 54. 35.) (2. 7. 12. 35. 43.)))return:'((4. 2. 6. 8. 2.) (0. 9. 6. 9. 6.) (0. 0. 1. 2. 5.) (0. 0. 0. 25. 15.) (0. 0. 0. 0. 4.))must return:'((2. 1. 3. 4. 1.) (0. 3. 2. 3. 2.) (0. 0. 1. 2. 5.) (0. 0. 0.  5.  3.) (0. 0. 0. 0. 2.))`

#### Lee Mac

• Seagull
• Posts: 12943
• London, England
##### Re: Challenge: Cholesky decomposition
« Reply #7 on: March 10, 2010, 05:10:51 PM »
Hello Lee Mac
You turned interesting program, but there is a small mistake ...
For the test case:
Code: [Select]
`(LM_Cholesky '((4. 2. 6. 8. 2.) (2. 10. 9. 13. 7.) (6. 9. 14. 20. 12.) (8. 13. 20. 54. 35.) (2. 7. 12. 35. 43.)))return:'((4. 2. 6. 8. 2.) (0. 9. 6. 9. 6.) (0. 0. 1. 2. 5.) (0. 0. 0. 25. 15.) (0. 0. 0. 0. 4.))must return:'((2. 1. 3. 4. 1.) (0. 3. 2. 3. 2.) (0. 0. 1. 2. 5.) (0. 0. 0.  5.  3.) (0. 0. 0. 0. 2.))`

Ah yes, I'm out by a factor. I think I actually did LU decomposition instead of Cholesky too...

#### Lee Mac

• Seagull
• Posts: 12943
• London, England
##### Re: Challenge: Cholesky decomposition
« Reply #8 on: March 10, 2010, 06:19:11 PM »
Slightly better, but still not what you want...

Code: [Select]
`(defun LM_Cholesky (mat / v+v sxv put i j)  ;; Lee Mac  ~  10.03.10  (defun v-v (v1 v2)    (mapcar (function -) v1 v2))  (defun sxv (s v1)    (mapcar (function (lambda (x) (* s x))) v1))  (defun put (x n lst / k)    (setq k -1)    (mapcar (function (lambda (y) (if (= (setq k (1+ k)) n) x y))) lst))      (setq i -1)     (mapcar    (function      (lambda (row / r)        (setq i (1+ i) j -1)        (repeat i                    (setq mat                 (put                      (setq row                             (v-v row                                     (setq r                                            (sxv (/ (nth (setq j (1+ j)) row)                                                    (float (nth j (nth j mat))))                                                 (nth j mat)))))                         i mat)          )          (setq mat (put r j mat))        )      )    )    mat)     mat)`

#### ElpanovEvgeniy

• Water Moccasin
• Posts: 1569
• Moscow (Russia)
##### Re: Challenge: Cholesky decomposition
« Reply #9 on: March 11, 2010, 12:32:04 AM »
return:
Code: [Select]
`'((4. 2. 6. 8. 2.) (0. 9. 6. 9. 6.) (0. 0. 1. 2. 5.) (0. 0. 0. 25. 15.) (0. 0. 0. 0. 4.))`
you just forgot:
Code: [Select]
`(append (list (setq i (sqrt 4.))) (mapcar '(lambda (a) (/ a i)) '(2. 6. 8. 2.)))(append (list 0.(setq i (sqrt 9.))) (mapcar '(lambda (a) (/ a i)) '( 6. 9. 6.)))(append (list 0. 0.(setq i (sqrt 1.))) (mapcar '(lambda (a) (/ a i)) '(2. 5.)))(append (list 0. 0. 0.(setq i (sqrt 25.))) (mapcar '(lambda (a) (/ a i)) '(15.)))(append (list 0. 0. 0. 0.(setq i (sqrt 4.))) (mapcar '(lambda (a) (/ a i)) '()))`

#### Lee Mac

• Seagull
• Posts: 12943
• London, England
##### Re: Challenge: Cholesky decomposition
« Reply #10 on: March 11, 2010, 06:23:40 AM »
Ahh! Thanks for that Evgeniy

Code: [Select]
`(defun LM_Cholesky (mat / v+v sxv put i j)  ;; Lee Mac  ~  10.03.10  (defun v-v (v1 v2)    (mapcar (function -) v1 v2))  (defun sxv (s v1)    (mapcar (function (lambda (x) (* s x))) v1))  (defun put (x n lst / k)    (setq k -1)    (mapcar (function (lambda (y) (if (= (setq k (1+ k)) n) x y))) lst))      (setq i -1)     (mapcar    (function      (lambda (row / r)        (setq i (1+ i) j -1)        (repeat i                    (setq mat            (put              (setq row                (v-v row                  (sxv (/ (nth (setq j (1+ j)) row)                          (float (nth j (nth j mat))))                                              (nth j mat)                  )                )              )            i mat)          )        )      )    )    mat)  (setq i -1)  (mapcar    (function      (lambda (row / rt)        (setq rt (sqrt (nth (setq i (1+ i)) row)))        (mapcar          (function            (lambda (element) (/ element rt)))          row)))    mat))`

#### ElpanovEvgeniy

• Water Moccasin
• Posts: 1569
• Moscow (Russia)
##### Re: Challenge: Cholesky decomposition
« Reply #11 on: March 11, 2010, 06:57:52 AM »
Ahh! Thanks for that Evgeniy

now great!

Very strange, the topic is not popular ...

#### Lee Mac

• Seagull
• Posts: 12943
• London, England
##### Re: Challenge: Cholesky decomposition
« Reply #12 on: March 11, 2010, 08:29:42 AM »
now great!

Thanks

Quote from: ElpanovEvgeniy
Very strange, the topic is not popular ...

Well it is quite a pure math topic - quite specialist

#### ElpanovEvgeniy

• Water Moccasin
• Posts: 1569
• Moscow (Russia)
##### Re: Challenge: Cholesky decomposition
« Reply #13 on: March 11, 2010, 12:37:40 PM »
Well it is quite a pure math topic - quite specialist

my version:
Code: [Select]
`(defun ai (a b) (sqrt (apply (function -) (cons a (mapcar (function *) b b)))))(defun aii (a b c) (cons a       (mapcar        (function         (lambda (b d) (/ (apply (function -) (cons b (mapcar (function *) (car c) d))) a))        )        b        (cdr c)       ) ))(defun cholesky (a l / ll) (cond  ((not l) (cholesky a (list (cons 0. (mapcar (function (lambda (a) 0.)) (car a))))))  (a   (cons    (setq ll (aii (ai (caar a) (mapcar (function car) l))                  (cdar a)                  (apply (function mapcar) (cons 'list l))             )    )    (mapcar (function (lambda (a) (cons 0. a)))            (cholesky (mapcar (function cdr) (cdr a)) (mapcar (function cdr) (cons ll l)))    )   )  )  (nil) ))`
Test:
Code: [Select]
`(cholesky matr nil)`
« Last Edit: March 11, 2010, 12:41:21 PM by ElpanovEvgeniy »

#### Lee Mac

• Seagull
• Posts: 12943
• London, England
##### Re: Challenge: Cholesky decomposition
« Reply #14 on: March 11, 2010, 12:51:17 PM »
Absolutely brilliant Evgeniy

#### pkohut

• Guest
##### Re: Challenge: Cholesky decomposition
« Reply #15 on: March 11, 2010, 12:55:03 PM »
Very strange, the topic is not popular ...

This challenge went beyond what I know, and the code from the links Swift provided is highly optimized.  So it's watch and learn time

#### T.Willey

• Needs a day job
• Posts: 5251
##### Re: Challenge: Cholesky decomposition
« Reply #16 on: March 11, 2010, 12:57:41 PM »
Very strange, the topic is not popular ...

This challenge went beyond what I know, and the code from the links Swift provided is highly optimized.  So it's watch and learn time
x2
Tim

I don't want to ' end-up ', I want to ' become '. - Me

#### ElpanovEvgeniy

• Water Moccasin
• Posts: 1569
• Moscow (Russia)
##### Re: Challenge: Cholesky decomposition
« Reply #17 on: March 11, 2010, 01:00:06 PM »
This code can be useful if you use the method of least squares
http://en.wikipedia.org/wiki/Least_squares

#### ElpanovEvgeniy

• Water Moccasin
• Posts: 1569
• Moscow (Russia)
##### Re: Challenge: Cholesky decomposition
« Reply #18 on: March 11, 2010, 01:01:47 PM »
Very strange, the topic is not popular ...

This challenge went beyond what I know, and the code from the links Swift provided is highly optimized.  So it's watch and learn time
x2

I was too early to show my version?

#### T.Willey

• Needs a day job
• Posts: 5251
##### Re: Challenge: Cholesky decomposition
« Reply #19 on: March 11, 2010, 01:17:31 PM »
Very strange, the topic is not popular ...

This challenge went beyond what I know, and the code from the links Swift provided is highly optimized.  So it's watch and learn time
x2

I was too early to show my version?

Not for me, as I don't really have time to study matrix math right now.
Tim

I don't want to ' end-up ', I want to ' become '. - Me

#### Swift

• Swamp Rat
• Posts: 596
##### Re: Challenge: Cholesky decomposition
« Reply #20 on: March 11, 2010, 01:39:06 PM »
This code can be useful if you use the method of least squares
http://en.wikipedia.org/wiki/Least_squares

Depending on your purpose you may want to look at LU Decomposition. Cornbread and I did a lot of work with this several years ago for some least squares coordinate transformations and I am wanting to think that is was easier to calculate the inverse of a square matrix with LU Decomposition.

#### Lee Mac

• Seagull
• Posts: 12943
• London, England
##### Re: Challenge: Cholesky decomposition
« Reply #21 on: March 11, 2010, 01:46:08 PM »
This code can be useful if you use the method of least squares
http://en.wikipedia.org/wiki/Least_squares

Depending on your purpose you may want to look at LU Decomposition. Cornbread and I did a lot of work with this several years ago for some least squares coordinate transformations and I am wanting to think that is was easier to calculate the inverse of a square matrix with LU Decomposition.

That's the method I used     But, according to Evgeniy's link, the Cholesky algorithm is meant to be much more efficient than LU decomposition...

#### Swift

• Swamp Rat
• Posts: 596
##### Re: Challenge: Cholesky decomposition
« Reply #22 on: March 11, 2010, 01:56:01 PM »
That's the method I used     But, according to Evgeniy's link, the Cholesky algorithm is meant to be much more efficient than LU decomposition...

Go back and look at the Wiki page
Quote
Cholesky decomposition or Cholesky triangle is a decomposition of a symmetric, positive-definite matrix into the product of a lower triangular matrix and its conjugate transpose

For the types of least squares solutions I work, you almost never have a symmetrical matrix where A = A^T

LU Decomposition doesn't have this limitation.

#### ElpanovEvgeniy

• Water Moccasin
• Posts: 1569
• Moscow (Russia)
##### Re: Challenge: Cholesky decomposition
« Reply #23 on: March 11, 2010, 01:56:24 PM »
Depending on your purpose you may want to look at LU Decomposition. Cornbread and I did a lot of work with this several years ago for some least squares coordinate transformations and I am wanting to think that is was easier to calculate the inverse of a square matrix with LU Decomposition.

I am well acquainted with the LU Decomposition.

#### Swift

• Swamp Rat
• Posts: 596
##### Re: Challenge: Cholesky decomposition
« Reply #24 on: March 11, 2010, 02:00:59 PM »
Very strange, the topic is not popular ...

This challenge went beyond what I know, and the code from the links Swift provided is highly optimized.  So it's watch and learn time

The Numerical Recipes book is a great book, full of lots of information, but be sure to read the copyrights before you buy it or the source code. Where as the TNT is published by the National Institute of Standards and Technology and is government funded and "public domain" as I understand it.

#### ElpanovEvgeniy

• Water Moccasin
• Posts: 1569
• Moscow (Russia)
##### Re: Challenge: Cholesky decomposition
« Reply #25 on: March 11, 2010, 02:11:29 PM »
For the types of least squares solutions I work, you almost never have a symmetrical matrix where A = A^T

Such matrices often arise when using the method of least squares

LU Decomposition doesn't have this limitation.

I agree with you, but it is more convenient to have two tools, but not one ...

#### Swift

• Swamp Rat
• Posts: 596
##### Re: Challenge: Cholesky decomposition
« Reply #26 on: March 11, 2010, 02:14:21 PM »
Depending on your purpose you may want to look at LU Decomposition. Cornbread and I did a lot of work with this several years ago for some least squares coordinate transformations and I am wanting to think that is was easier to calculate the inverse of a square matrix with LU Decomposition.

I am well acquainted with the LU Decomposition.

So it seems you are. I just made that comment to perhaps save you time. I spent a lot of time chasing my tail trying to develop a matrix inversion routine before I found the LU Decomposition.

Agreed on having two tools.

#### LE3

• Guest
##### Re: Challenge: Cholesky decomposition
« Reply #27 on: March 11, 2010, 02:49:40 PM »
from a non math persona...
and nothing to offer, but I have implemented the Gaussian elimination at least for two of my apps, inside of gbpoly to do the pivoting onto the recognition of closed circuits and the obb extraction both in C++.

and must say that do not recall hearing the cholesky named from my school days.... but remember something about 'factorizacion' (factorization) but that's all.
« Last Edit: March 11, 2010, 02:54:10 PM by LE3 »