least-squares estimate
Evan Monroig <[email protected]> Tue, 22 Apr 2008 14:43:15 +0900
| Newsgroups | gmane.lisp.matlisp.user |
|---|---|
| Message-ID | <[email protected]> |
--=-=-= Hi, I found the need to find the least-squares solution of a classical equation Ax = b, so I added a binding for the lapack function DGELSY. This is the one that solves rectangular real linear least squares problems using a complete orthogonal factorization (see [1] and [2]). Please find a patch attached. I authorize you to incorporate it to matlisp as is or after modification, using the same license as other parts of matlisp. Cheers, Evan [1] http://www.netlib.org/lapack/lug/node27.html [2] http://www.netlib.org/lapack/lug/node43.html --=-=-= Content-Type: text/x-diff Content-Disposition: attachment; filename=0001-packages.lisp-LAPACK-export-DGELSY.patch Content-Description: 0001-packages.lisp-LAPACK-export-DGELSY.patch >From 68827166d3fcd0f4cb58f51e40c87792617e7715 Mon Sep 17 00:00:00 2001 From: Evan Monroig <[email protected]> Date: Tue, 22 Apr 2008 14:29:47 +0900 Subject: [PATCH] * packages.lisp ("LAPACK"): export DGELSY ("MATLISP"): export GELSY! and GELSY * system.dcl (matlisp): added file src/gels.lisp * src/gels.lisp (gelsy!, gelsy): wrappers over DGELSY * matlisp.mk.in (LAPACK_OBJS): add dgelsy.o and dependencies * src/lapack.lisp (dgelsy): added - minimum-norm solution to a real linear least squares problem using a complete orthogonal factorization --- matlisp.mk.in | 5 ++- packages.lisp | 4 +- src/gels.lisp | 145 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/lapack.lisp | 139 ++++++++++++++++++++++++++++++++++++++++++++++++++++ system.dcl | 3 +- 5 files changed, 293 insertions(+), 3 deletions(-) create mode 100644 src/gels.lisp diff --git a/matlisp.mk.in b/matlisp.mk.in index b709036..3df9b96 100644 --- a/matlisp.mk.in +++ b/matlisp.mk.in @@ -139,7 +139,10 @@ LAPACK_OBJS = dlasq5.o dlasq6.o ieeeck.o zdrot.o dlabad.o \ dggbal.o dlansy.o dorg2r.o izmax1.o zgetf2.o zlasr.o zunmlq.o \ dgghrd.o dlanv2.o dorgbr.o zlassq.o zunmqr.o \ dhgeqz.o dlapy2.o dorghr.o \ - dhseqr.o dlapy3.o dorgl2.o zheev.o zlatrd.o $(NO_ATLAS_LAPACK_OBJS) + dhseqr.o dlapy3.o dorgl2.o zheev.o zlatrd.o \ + dlaic1.o dlatrz.o dlarz.o dlarzb.o dlarzt.o dgelsy.o dtzrzf.o \ + dormr3.o dormrz.o\ + $(NO_ATLAS_LAPACK_OBJS) # The FFT package DFFTPACK_OBJS = zfftb.o \ diff --git a/packages.lisp b/packages.lisp index 05b9eaa..7c6f343 100644 --- a/packages.lisp +++ b/packages.lisp @@ -161,7 +161,7 @@ "DGESV" "DGEEV" "DGETRF" "DGETRS" "DGESVD" "ZGESV" "ZGEEV" "ZGETRF" "ZGETRS" "ZGESVD" "DGEQRF" "ZGEQRF" "DGEQP3" "ZGEQP3" - "DORGQR" "ZUNGQR") + "DORGQR" "ZUNGQR" "DGELSY") (:documentation "LAPACK routines")) (defpackage "DFFTPACK" @@ -275,6 +275,8 @@ "FORTRAN-COMPLEX-MATRIX-INDEXING" "FORTRAN-MATRIX-INDEXING" "GEEV" + "GELSY!" + "GELSY" "GEMM!" "GEMM" "GESV!" diff --git a/src/gels.lisp b/src/gels.lisp new file mode 100644 index 0000000..0d8008d --- /dev/null +++ b/src/gels.lisp @@ -0,0 +1,145 @@ +;;; -*- Mode: lisp; Syntax: ansi-common-lisp; Package: :matlisp; Base: 10 -*- + +(in-package "MATLISP") + +(defgeneric gelsy! (a b rcond) + (:documentation "Destructive version of GELSY. See GELSY.")) + +(defgeneric gelsy (a b rcond) + (:documentation + " + Syntax + ======= + + (GELSY A B &key TOL) + + INPUT + ----- + A A Matlisp matrix of size M x N + B A Matlisp matrix of size M x P + RCOND A condition number + + OUTPUT + ------ + X A Matlisp matrix of size N x NRHS + RANK An integer + + Purpose + ======= + + Compute the minimum-norm solution to a real linear least + squares problem: + minimize || A * X - B || + using a complete orthogonal factorization of A. A is an M-by-N + matrix which may be rank-deficient. + + Several right hand side vectors b and solution vectors x can be + handled in a single call; they are stored as the columns of the + M-by-NRHS right hand side matrix B and the N-by-NRHS solution + matrix X. + + The routine first computes a QR factorization with column pivoting: + A * P = Q * [ R11 R12 ] + [ 0 R22 ] + with R11 defined as the largest leading submatrix whose estimated + condition number is less than 1/RCOND. The order of R11, RANK, + is the effective rank of A. + + Then, R22 is considered to be negligible, and R12 is annihilated + by orthogonal transformations from the right, arriving at the + complete orthogonal factorization: + A * P = Q * [ T11 0 ] * Z + [ 0 0 ] + The minimum-norm solution is then + X = P * Z' [ inv(T11)*Q1'*B ] + [ 0 ] + where Q1 consists of the first RANK columns of Q. + + This routine is basically identical to the original xGELSX except + three differences: + o The call to the subroutine xGEQPF has been substituted by the + the call to the subroutine xGEQP3. This subroutine is a Blas-3 + version of the QR factorization with column pivoting. + o Matrix B (the right hand side) is updated with Blas-3. + o The permutation of matrix B (the right hand side) is faster and + more simple. + + Further Details + =============== + + Based on contributions by + A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA + E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain + G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain + + ===================================================================== +")) + +(defun check-info (info function-name) + (unless (= info 0) + (error "~a: error in argument ~d" function-name (- info)))) + +(defun dgelsy-workspace-inquiry (m n nrhs a lda b ldb jpvt rcond rank) + (let ((work (allocate-real-store 1))) + (multiple-value-bind + (store-a store-b store-jpvt rank store-work info) + (lapack::dgelsy m n nrhs (store a) lda (store b) ldb jpvt rcond rank + work -1 0) + (declare (ignore store-a store-b store-jpvt rank store-work)) + (check-info info "dgelsy")) + (values (ceiling (realpart (aref work 0)))))) + +(defmethod gelsy! ((a real-matrix) (b real-matrix) rcond) + (let* ((m (nrows a)) + (n (ncols a)) + (nrhs (ncols b)) + (jpvt (allocate-integer4-store n 0)) + (lda m) + (ldb (max n m)) + (b-arg b)) + (when (and (< m n)) + ;; In this case we need to extend the matrix which stores B + ;; since it will be used to store the computation result + (setq b-arg (make-real-matrix n nrhs)) + (dotimes (i m) + (dotimes (j nrhs) + (setf (matrix-ref b-arg i j) + (matrix-ref b i j))))) + (let* ((lwork (dgelsy-workspace-inquiry m n nrhs a lda b-arg ldb jpvt + rcond 0)) + (work (allocate-real-store lwork))) + (assert (= m (nrows b))) + (multiple-value-bind + (store-a store-b store-jpvt rank store-work info) + (lapack::dgelsy m n nrhs (store a) lda (store b-arg) ldb jpvt rcond 0 + work lwork 0) + (declare (ignore store-a store-jpvt store-work)) + (check-info info "dgelsy") + (let ((x (make-real-matrix n nrhs))) + ;; extract the matrix X from B + (dotimes (i n) + (dotimes (j nrhs) + (setf (matrix-ref x i j) + (aref store-b (fortran-matrix-indexing i j n))))) + (values x rank)))))) + +(defmethod gelsy ((a real-matrix) (b real-matrix) rcond) + (gelsy! (copy a) (copy b) rcond)) + +;; Example + +#| + +(let* ((m 100) + (n 100) + (a (rand m n)) + (x (rand n 1)) + (b (m* a x)) + (eps (coerce (expt 2 -52) 'double-float)) + (rcond (* eps (max m n)))) + (multiple-value-bind (r1 rank) + (matlisp::gelsy a b rcond) + (list rank + (norm (m- x r1))))) + +|# diff --git a/src/lapack.lisp b/src/lapack.lisp index b49adff..b355c98 100644 --- a/src/lapack.lisp +++ b/src/lapack.lisp @@ -1510,3 +1510,142 @@ (lwork :integer :input) (rwork (* :double-float) :workspace-output) (info :integer :output)) + +(def-fortran-routine dgelsy :void + " + -- LAPACK driver routine (version 3.0) -- + Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., + Courant Institute, Argonne National Lab, and Rice University + June 30, 1999 + + Purpose + ======= + + DGELSY computes the minimum-norm solution to a real linear least + squares problem: + minimize || A * X - B || + using a complete orthogonal factorization of A. A is an M-by-N + matrix which may be rank-deficient. + + Several right hand side vectors b and solution vectors x can be + handled in a single call; they are stored as the columns of the + M-by-NRHS right hand side matrix B and the N-by-NRHS solution + matrix X. + + The routine first computes a QR factorization with column pivoting: + A * P = Q * [ R11 R12 ] + [ 0 R22 ] + with R11 defined as the largest leading submatrix whose estimated + condition number is less than 1/RCOND. The order of R11, RANK, + is the effective rank of A. + + Then, R22 is considered to be negligible, and R12 is annihilated + by orthogonal transformations from the right, arriving at the + complete orthogonal factorization: + A * P = Q * [ T11 0 ] * Z + [ 0 0 ] + The minimum-norm solution is then + X = P * Z' [ inv(T11)*Q1'*B ] + [ 0 ] + where Q1 consists of the first RANK columns of Q. + + This routine is basically identical to the original xGELSX except + three differences: + o The call to the subroutine xGEQPF has been substituted by the + the call to the subroutine xGEQP3. This subroutine is a Blas-3 + version of the QR factorization with column pivoting. + o Matrix B (the right hand side) is updated with Blas-3. + o The permutation of matrix B (the right hand side) is faster and + more simple. + + Arguments + ========= + + M (input) INTEGER + The number of rows of the matrix A. M >= 0. + + N (input) INTEGER + The number of columns of the matrix A. N >= 0. + + NRHS (input) INTEGER + The number of right hand sides, i.e., the number of + columns of matrices B and X. NRHS >= 0. + + A (input/output) DOUBLE PRECISION array, dimension (LDA,N) + On entry, the M-by-N matrix A. + On exit, A has been overwritten by details of its + complete orthogonal factorization. + + LDA (input) INTEGER + The leading dimension of the array A. LDA >= max(1,M). + + B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) + On entry, the M-by-NRHS right hand side matrix B. + On exit, the N-by-NRHS solution matrix X. + + LDB (input) INTEGER + The leading dimension of the array B. LDB >= max(1,M,N). + + JPVT (input/output) INTEGER array, dimension (N) + On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted + to the front of AP, otherwise column i is a free column. + On exit, if JPVT(i) = k, then the i-th column of AP + was the k-th column of A. + + RCOND (input) DOUBLE PRECISION + RCOND is used to determine the effective rank of A, which + is defined as the order of the largest leading triangular + submatrix R11 in the QR factorization with pivoting of A, + whose estimated condition number < 1/RCOND. + + RANK (output) INTEGER + The effective rank of A, i.e., the order of the submatrix + R11. This is the same as the order of the submatrix T11 + in the complete orthogonal factorization of A. + + WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) + On exit, if INFO = 0, WORK(1) returns the optimal LWORK. + + LWORK (input) INTEGER + The dimension of the array WORK. + The unblocked strategy requires that: + LWORK >= MAX( MN+3*N+1, 2*MN+NRHS ), + where MN = min( M, N ). + The block algorithm requires that: + LWORK >= MAX( MN+2*N+NB*(N+1), 2*MN+NB*NRHS ), + where NB is an upper bound on the blocksize returned + by ILAENV for the routines DGEQP3, DTZRZF, STZRQF, DORMQR, + and DORMRZ. + + If LWORK = -1, then a workspace query is assumed; the routine + only calculates the optimal size of the WORK array, returns + this value as the first entry of the WORK array, and no error + message related to LWORK is issued by XERBLA. + + INFO (output) INTEGER + = 0: successful exit + < 0: If INFO = -i, the i-th argument had an illegal value. + + Further Details + =============== + + Based on contributions by + A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA + E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain + G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain + + ===================================================================== +" + (m :integer :input) + (n :integer :input) + (nrhs :integer :input) + (a (* :double-float) :input-output) + (lda :integer :input) + (b (* :double-float) :input-output) + (ldb :integer :input) + (jpvt (* :integer) :input-output) + (rcond :double-float :input) + (rank :integer :output) + (work (* :double-float) :workspace-output) + (lwork :integer :input) + (info :integer :output)) diff --git a/system.dcl b/system.dcl index 329020c..de1c231 100644 --- a/system.dcl +++ b/system.dcl @@ -205,7 +205,8 @@ :depends-on ("foreign-interface" "foreign-functions" "matlisp-essentials") - :components ("gesv" + :components ("gels" + "gesv" "geev" "getrf" "getrs")) -- 1.5.2.5 --=-=-= Content-Type: text/plain; charset="us-ascii" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit Content-Disposition: inline ------------------------------------------------------------------------- This SF.net email is sponsored by the 2008 JavaOne(SM) Conference Don't miss this year's exciting event. There's still time to save $100. Use priority code J8TL2D2. http://ad.doubleclick.net/clk;198757673;13503038;p?http://java.sun.com/javaone --=-=-= Content-Type: text/plain; charset="us-ascii" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit Content-Disposition: inline _______________________________________________ Matlisp-users mailing list [email protected] https://lists.sourceforge.net/lists/listinfo/matlisp-users --=-=-=--