Commit 03ad93e0 authored by Liam Healy's avatar Liam Healy
Browse files

SVD: More optional arguments, add tests

For the singular value decomposition functions, make the S vector and
V matrix optional arguments, as they are quantities returned by the
functions.  Add tests translated from the GSL tests.
parent 2965e775
;; Definition of GSLL system
;; Liam Healy
;; Time-stamp: <2009-08-26 21:04:08EDT gsll-tests.asd>
;; Time-stamp: <2009-09-18 16:26:18EDT gsll-tests.asd>
(asdf:defsystem "gsll-tests"
:name "gsll-tests"
......@@ -159,6 +159,7 @@
(:file "sort-vector-smallest-index")
(:file "sort-vector-smallest")
(:file "spherical-vector")
(:file "svd")
(:file "swap-columns")
(:file "swap-elements")
(:file "swap-row-column")
......
;; Singular Value Decomposition
;; Liam Healy, Tue May 2 2006 - 12:15
;; Time-stamp: <2009-02-23 20:52:59EST svd.lisp>
;; Time-stamp: <2009-09-18 16:07:28EDT svd.lisp>
;; $Id$
(in-package :gsl)
;;; /usr/include/gsl/gsl_linalg.h
;;; FDL
;;; A general rectangular M-by-N matrix A has a
;;; singular value decomposition (svd) into the product of an
;;; M-by-N orthogonal matrix U, an N-by- diagonal matrix of
;;; M-by-N orthogonal matrix U, an N-by-N diagonal matrix of
;;; singular values S and the transpose of an
;;; N-by-N orthogonal square matrix V, A = U S V^T
;;; The singular values sigma_i = S_{ii} are all non-negative and are
......@@ -26,7 +28,10 @@
;;; tolerance.
(defmfun SV-decomposition
(A S V &optional (work (make-marray 'double-float :dimensions (dim1 A))))
(A &optional
(S (make-marray 'double-float :dimensions (dim1 A)))
(V (make-marray 'double-float :dimensions (list (dim1 A) (dim1 A))))
(work (make-marray 'double-float :dimensions (dim1 A))))
"gsl_linalg_SV_decomp"
(((mpointer A) :pointer) ((mpointer V) :pointer)
((mpointer S) :pointer) ((mpointer work) :pointer))
......@@ -46,9 +51,12 @@
This routine uses the Golub-Reinsch SVD algorithm.")
(defmfun SV-modified-decomposition
(A S V
&optional (X (make-marray 'double-float :dimensions (list (dim1 A) (dim1 A))))
(work (make-marray 'double-float :dimensions (dim1 A))))
(A
&optional
(S (make-marray 'double-float :dimensions (dim1 A)))
(V (make-marray 'double-float :dimensions (list (dim1 A) (dim1 A))))
(X (make-marray 'double-float :dimensions (list (dim1 A) (dim1 A))))
(work (make-marray 'double-float :dimensions (dim1 A))))
"gsl_linalg_SV_decomp_mod"
(((mpointer A) :pointer) ((mpointer X) :pointer)
((mpointer V) :pointer)
......@@ -61,7 +69,10 @@
faster for M >> N. It requires the vector work of length N and the
N-by-N matrix X as additional working space.")
(defmfun SV-jacobi-decomposition (A S V)
(defmfun SV-jacobi-decomposition
(A &optional
(S (make-marray 'double-float :dimensions (dim1 A)))
(V (make-marray 'double-float :dimensions (list (dim1 A) (dim1 A)))))
"gsl_linalg_SV_decomp_jacobi"
(((mpointer A) :pointer) ((mpointer V) :pointer)
((mpointer S) :pointer))
......@@ -95,3 +106,44 @@
In the over-determined case where A has more rows than columns the
system is solved in the least squares sense, returning the solution
x which minimizes ||A x - b||_2.")
;;; Examples and unit test, from linalg/test.c
;;; These are general to all the linear solver techniques, so
;;; more tests need to be made.
(defun create-hilbert-matrix (dim)
"Make Hilbert matrix used to test linear algebra functions."
(let ((matrix (make-marray 'double-float :dimensions (list dim dim))))
(dotimes (i dim matrix)
(dotimes (j dim)
(setf (maref matrix i j) (coerce (/ (+ 1 i j)) 'double-float))))))
(defun create-vandermonde-matrix (dim)
"Make Van der Monde matrix used to test linear algebra functions."
(let ((matrix (make-marray 'double-float :dimensions (list dim dim))))
(dotimes (i dim matrix)
(dotimes (j dim)
(setf (maref matrix i j)
(coerce (expt (1+ i) (- dim j 1)) 'double-float))))))
(defun test-sv-solve-dim (matrix)
"Solve the linear equation using SVD with the supplied matrix and
a right-hand side vector which is the reciprocal of one more than
the index."
(let* ((dim (dim0 matrix))
(rhs (make-marray 'double-float :dimensions dim)))
(dotimes (i dim)
(setf (maref rhs i) (coerce (1+ i) 'double-float)))
(multiple-value-bind (u q d)
(SV-decomposition (copy matrix))
(SV-solve u q d rhs))))
(save-test svd
(test-sv-solve-dim (create-hilbert-matrix 2))
(test-sv-solve-dim (create-hilbert-matrix 3))
(test-sv-solve-dim (create-hilbert-matrix 4))
(test-sv-solve-dim (create-hilbert-matrix 12))
(test-sv-solve-dim (create-vandermonde-matrix 2))
(test-sv-solve-dim (create-vandermonde-matrix 3))
(test-sv-solve-dim (create-vandermonde-matrix 4))
(test-sv-solve-dim (create-vandermonde-matrix 12)))
;; Regression test SVD for GSLL, aputomatically generated
(in-package :gsl)
;;; Answers inserted from linalg/test.c
;;; GSL has #define GSL_DBL_EPSILON 2.2204460492503131e-16
;;; which is 2x what double-float-epsilon is.
(LISP-UNIT:DEFINE-TEST SVD
(let ((lisp-unit:*epsilon* (* 2 16 double-float-epsilon)))
(LISP-UNIT:ASSERT-NUMERICAL-EQUAL
(LIST
(MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS '(-8.0d0 18.0d0)))
(MULTIPLE-VALUE-LIST
(TEST-SV-SOLVE-DIM (CREATE-HILBERT-MATRIX 2)))))
(let ((lisp-unit:*epsilon* (* 2 128 double-float-epsilon)))
(LISP-UNIT:ASSERT-NUMERICAL-EQUAL
(LIST
(MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS '(27.0d0 -192.0d0 210.0d0)))
(MULTIPLE-VALUE-LIST
(TEST-SV-SOLVE-DIM (CREATE-HILBERT-MATRIX 3)))))
(let ((lisp-unit:*epsilon* (* 2 2048 double-float-epsilon)))
(LISP-UNIT:ASSERT-NUMERICAL-EQUAL
(LIST
(MAKE-MARRAY 'DOUBLE-FLOAT
:INITIAL-CONTENTS '(-64.0d0 900.0d0 -2520.0d0 1820.0d0)))
(MULTIPLE-VALUE-LIST
(TEST-SV-SOLVE-DIM (CREATE-HILBERT-MATRIX 4)))))
(let ((lisp-unit:*epsilon* 0.5d0))
(LISP-UNIT:ASSERT-NUMERICAL-EQUAL
(LIST
(MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS
'(-1728.0d0 245388.0d0 -8528520.0d0
127026900.0d0 -1009008000.0d0 4768571808.0d0
-14202796608.0d0 27336497760.0d0 -33921201600.0d0
26189163000.0d0 -11437874448.0d0 2157916488.0d0)))
(MULTIPLE-VALUE-LIST
(TEST-SV-SOLVE-DIM (CREATE-HILBERT-MATRIX 12)))))
(let ((lisp-unit:*epsilon* (* 2 64 double-float-epsilon)))
(LISP-UNIT:ASSERT-NUMERICAL-EQUAL
(LIST (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS '(1.0d0 0.0d0)))
(MULTIPLE-VALUE-LIST
(TEST-SV-SOLVE-DIM (CREATE-VANDERMONDE-MATRIX 2)))))
(let ((lisp-unit:*epsilon* (* 2 64 double-float-epsilon)))
(LISP-UNIT:ASSERT-NUMERICAL-EQUAL
(LIST (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS '(0.0d0 1.0d0 0.0d0)))
(MULTIPLE-VALUE-LIST
(TEST-SV-SOLVE-DIM (CREATE-VANDERMONDE-MATRIX 3)))))
(let ((lisp-unit:*epsilon* (* 2 1024 double-float-epsilon)))
(LISP-UNIT:ASSERT-NUMERICAL-EQUAL
(LIST (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS '(0.0d0 0.0d0 1.0d0 0.0d0)))
(MULTIPLE-VALUE-LIST
(TEST-SV-SOLVE-DIM (CREATE-VANDERMONDE-MATRIX 4)))))
(let ((lisp-unit:*epsilon* 0.05d0))
(LISP-UNIT:ASSERT-NUMERICAL-EQUAL
(LIST
(MAKE-MARRAY
'DOUBLE-FLOAT :INITIAL-CONTENTS
'(0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0 0.0d0)))
(MULTIPLE-VALUE-LIST
(TEST-SV-SOLVE-DIM (CREATE-VANDERMONDE-MATRIX 12))))))
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment