cholesky.lisp 4.62 KB
Newer Older
liam's avatar
liam committed
1 2
;; Cholesky Decomposition
;; Liam Healy, Wed May  3 2006 - 16:38
3
;; Time-stamp: <2011-02-19 13:08:24EST cholesky.lisp>
4
;;
5
;; Copyright 2006, 2007, 2008, 2009, 2011 Liam M. Healy
6 7 8 9 10 11 12 13 14 15 16 17 18 19
;; Distributed under the terms of the GNU General Public License
;;
;; This program is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this program.  If not, see <http://www.gnu.org/licenses/>.
liam's avatar
liam committed
20 21 22

(in-package :gsl)

liam's avatar
liam committed
23 24 25 26
;;; FDL
;;; A symmetric, positive definite square matrix A has a Cholesky
;;; decomposition into a product of a lower triangular matrix L and
;;; its transpose L^T,
liam's avatar
liam committed
27 28 29 30
;;; A = L L^T
;;; This is sometimes referred to as taking the square-root of a matrix. The
;;; Cholesky decomposition can only be carried out when all the eigenvalues
;;; of the matrix are positive.  This decomposition can be used to convert
liam's avatar
liam committed
31 32
;;; the linear system A x = b into a pair of triangular systems
;;; (L y = b, L^T x = y), which can be solved by forward and
liam's avatar
liam committed
33 34
;;; back-substitution.

35 36
;;; GSL version 1.10 introduced functions for complex matrices.

37
(defmfun cholesky-decomposition ((A grid:matrix))
38 39 40 41
  ("gsl_linalg" :complex "_cholesky_decomp")
  (((mpointer A) :pointer))
  :definition :generic
  :element-types :doubles
42 43
  :inputs (A)
  :outputs (A)
liam's avatar
liam committed
44
  :documentation			; FDL
45 46 47 48 49 50 51
  "Factorize the positive-definite square matrix A into the Cholesky
  decomposition A = L L^T (real) or A = L L^H (complex).  On output the
  diagonal and lower triangular part of the input matrix A contain the
  matrix L.  The upper triangular part of the input matrix contains
  L^T, the diagonal terms being identical for both L and L^T.  If the
  matrix is not positive-definite then the decomposition will fail,
  returning the error input-domain.")
liam's avatar
liam committed
52

53
(defmfun cholesky-solve
54
    ((A grid:matrix) (b vector) &optional x-spec
55
     &aux
56
     (x (grid:ensure-foreign-array x-spec (grid:dimensions b) t)))
57 58 59 60 61 62 63
  (("gsl_linalg" :complex "_cholesky_svx")
   ("gsl_linalg" :complex "_cholesky_solve"))
  ((((mpointer A) :pointer) ((mpointer b) :pointer))
   (((mpointer A) :pointer) ((mpointer b) :pointer)
    ((mpointer x) :pointer)))
  :definition :generic
  :element-types :doubles
64 65
  :inputs (A b)
  :outputs (x)
66
  :return ((or x b))
liam's avatar
liam committed
67 68
  :documentation			; FDL
  "Solve the system A x = b using the Cholesky
69
  decomposition of A into the matrix given by
70 71
  #'cholesky-decomposition.  If x-spec is NIL (default), the solution
  will replace b.  If x-spec is T, then an array will be created and the
72
  solution returned in it.  If x-spec is a grid:foreign-array, the solution will
73
  be returned in it.")
74

Liam Healy's avatar
Liam Healy committed
75 76 77 78 79 80 81 82 83 84 85
(defmfun cholesky-invert (cholesky)
  "gsl_linalg_cholesky_invert"
  (((mpointer cholesky) :pointer))
  :inputs (cholesky)
  :outputs (cholesky)
  :gsl-version (1 12)
  :documentation
  "Compute the inverse of the matrix cholesky which must have been
  previously computed by #'cholesky-decomposition. The inverse of the
  original matrix is stored in cholesky on output.")

86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
;;; Examples and unit test, from linalg/test.c

(defun test-cholesky-solve-dim (matrix)
  "Solve the linear equation using Cholesky with the supplied matrix and
   a right-hand side vector which is the reciprocal of one more than
   the index."
  (cholesky-solve
   (cholesky-decomposition (copy matrix))
   (create-rhs-vector (dim0 matrix))))

(defun test-cholesky-decomp-dim (matrix)
  "Decompose using Cholesky and then multiply."
  (let ((decomp (cholesky-decomposition (copy matrix))))
    (dotimes (row (dim0 matrix) decomp)
      (loop for col from (1+ row) below (dim1 matrix) do
Liam Healy's avatar
Liam Healy committed
101
	   (setf (grid:gref decomp row col) 0.0d0)))
102 103
    (matrix-product decomp decomp nil 1.0d0 0.0d0 :notrans :trans)))

Liam Healy's avatar
Liam Healy committed
104 105 106 107
(defun test-cholesky-invert-dim (matrix)
  "Invert using Cholesky decomposition"
  (cholesky-invert (cholesky-decomposition (copy matrix))))

108 109 110 111 112 113 114 115
(save-test cholesky
 (test-cholesky-solve-dim *hilb2*)
 (test-cholesky-solve-dim *hilb3*)
 (test-cholesky-solve-dim *hilb4*)
 (test-cholesky-solve-dim *hilb12*)
 (test-cholesky-decomp-dim *hilb2*)
 (test-cholesky-decomp-dim *hilb3*)
 (test-cholesky-decomp-dim *hilb4*)
Liam Healy's avatar
Liam Healy committed
116 117 118 119 120
 (test-cholesky-decomp-dim *hilb12*)
 (test-cholesky-invert-dim *hilb2*)
 (test-cholesky-invert-dim *hilb3*)
 (test-cholesky-invert-dim *hilb4*)
 (test-cholesky-invert-dim *hilb12*))