Topic: Challenge: Cholesky decomposition

[ElpanovEvgeniy]

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]

Great challenge Evgeniy

[Swift]

I am suppose to be working on a pay application submittal but...

Look here
http://math.nist.gov/tnt/index.html
http://math.nist.gov/tnt/overview.html
http://math.nist.gov/tnt/jama_doxygen/class_JAMA__Cholesky.html

and here ( This is a great book ! )
http://www.nrbook.com/a/bookcpdf/c2-9.pdf

[pkohut]

Oh my, nice.

Can someone translate this for me? It's from the wiki link.

The algorithms described below all involve about n³/3 FLOPs, where n is the size of the matrix A

I don't have anything to reference it against.

[Swift]

Oh my, nice.

Can someone translate this for me? It's from the wiki link.

The algorithms described below all involve about n³/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]

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)

         


[ElpanovEvgeniy]

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]

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]

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]

your previous version:
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]

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]

Ahh! Thanks for that Evgeniy 

now great!

Very strange, the topic is not popular ...
I do not know already you can to spread its version or wait for additional variants?

[Lee Mac]

now great!

Thanks 

Very strange, the topic is not popular ...

Well it is quite a pure math topic - quite specialist 

[ElpanovEvgeniy]

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)

[Lee Mac]

Absolutely brilliant Evgeniy