cholesky.lisp 4.62 KB
 liam committed Feb 18, 2008 1 2 ``````;; Cholesky Decomposition ;; Liam Healy, Wed May 3 2006 - 16:38 `````` Liam M. Healy committed Feb 19, 2011 3 ``````;; Time-stamp: <2011-02-19 13:08:24EST cholesky.lisp> `````` Liam Healy committed Dec 27, 2009 4 ``````;; `````` Liam M. Healy committed Jan 10, 2011 5 ``````;; Copyright 2006, 2007, 2008, 2009, 2011 Liam M. Healy `````` Liam Healy committed Dec 27, 2009 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 . `````` liam committed May 03, 2006 20 21 22 `````` (in-package :gsl) `````` liam committed Feb 18, 2008 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 committed May 03, 2006 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 committed Feb 18, 2008 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 committed May 03, 2006 33 34 ``````;;; back-substitution. `````` Liam Healy committed Feb 16, 2009 35 36 ``````;;; GSL version 1.10 introduced functions for complex matrices. `````` Liam M. Healy committed Jan 12, 2011 37 ``````(defmfun cholesky-decomposition ((A grid:matrix)) `````` Liam Healy committed Feb 16, 2009 38 39 40 41 `````` ("gsl_linalg" :complex "_cholesky_decomp") (((mpointer A) :pointer)) :definition :generic :element-types :doubles `````` Liam Healy committed Aug 11, 2008 42 43 `````` :inputs (A) :outputs (A) `````` liam committed Feb 18, 2008 44 `````` :documentation ; FDL `````` Liam Healy committed Feb 16, 2009 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 committed May 03, 2006 52 `````` `````` Liam Healy committed Feb 16, 2009 53 ``````(defmfun cholesky-solve `````` Liam M. Healy committed Jan 12, 2011 54 `````` ((A grid:matrix) (b vector) &optional x-spec `````` Liam Healy committed Feb 17, 2009 55 `````` &aux `````` Liam M. Healy committed Feb 19, 2011 56 `````` (x (grid:ensure-foreign-array x-spec (grid:dimensions b) t))) `````` Liam Healy committed Feb 16, 2009 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 `````` Liam Healy committed Aug 11, 2008 64 65 `````` :inputs (A b) :outputs (x) `````` Liam Healy committed Feb 17, 2009 66 `````` :return ((or x b)) `````` liam committed Feb 18, 2008 67 68 `````` :documentation ; FDL "Solve the system A x = b using the Cholesky `````` Liam Healy committed Feb 16, 2009 69 `````` decomposition of A into the matrix given by `````` Liam Healy committed Feb 17, 2009 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 `````` Liam Healy committed Jun 29, 2010 72 `````` solution returned in it. If x-spec is a grid:foreign-array, the solution will `````` Liam Healy committed Feb 17, 2009 73 `````` be returned in it.") `````` Liam Healy committed Sep 26, 2009 74 `````` `````` Liam Healy committed Jan 18, 2010 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.") `````` Liam Healy committed Sep 26, 2009 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 committed Jun 30, 2010 101 `````` (setf (grid:gref decomp row col) 0.0d0))) `````` Liam Healy committed Sep 26, 2009 102 103 `````` (matrix-product decomp decomp nil 1.0d0 0.0d0 :notrans :trans))) `````` Liam Healy committed Jan 18, 2010 104 105 106 107 ``````(defun test-cholesky-invert-dim (matrix) "Invert using Cholesky decomposition" (cholesky-invert (cholesky-decomposition (copy matrix)))) `````` Liam Healy committed Sep 26, 2009 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 committed Jan 18, 2010 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*))``````