Cholesky matrix decomposition

cholesky-matrix-decomp.sps

;;; Based on the following post:
;;; http://groups.google.com/group/comp.lang.scheme/browse_thread/thread/78170599ca8be8f6#

#!r6rs
(import (rnrs))

(define (Cholesky:decomp P)
  (Cholesky:make-square
   (let iter ((i 0) (j 0) (L '()))
     (if (>= i (length P))
	 L
       (iter (if (= i j) (+ 1 i) i)
	     (if (= i j) 0 (+ 1 j))
	     (Cholesky:add-element P L i j))))))

(define (Cholesky:add-element P L i j)
  (if (zero? j)
      (append L `((,(Cholesky:make-element P L i j))))
    (append (Cholesky:smaller L)
	    (list (append
		   (Cholesky:last-row L)
		   (list (Cholesky:make-element P L i j)))))))

(define (Cholesky:make-element P L i j)
  (let ((x (- (matrix:element P i j)
	      (Cholesky:partial-sum L i j))))
    (if (= i j)
	(sqrt x)
      (/ x (matrix:element L j j)))))

(define (Cholesky:partial-sum L i j)
  (let loop ((k j))
    (case k
      ((0) 0)
      ((1) (* (matrix:element L i 0)
	      (matrix:element L j 0)))
      (else
       (+ (* (matrix:element L i k)
	     (matrix:element L j k))
	  (loop (- k 1)))))))

(define (matrix:element A i j)
  (list-ref (list-ref A i) j))

(define (Cholesky:last-row L)
  (car (reverse L)))

(define (Cholesky:smaller P)
  (if (null? (cdr P))
      '()
    (reverse (cdr (reverse P)))))

(define (zero-vector n)
  (if (zero? n)
      '()
    (cons 0 (zero-vector (- n 1)))))

(define (Cholesky:make-square L)
  (map (lambda (v)
	 (append v (zero-vector (- (length L) (length v)))))
    L))

;;; ----------

(define A '((2 -2)
	    (-2 5)))

(write (Cholesky:decomp A))
  ;; => '(( 1.414213562 0.0        )
  ;;      (-1.414213562 1.732050808))

;;; end of file