;; QR decomposition ;; Liam Healy 2008-02-17 11:05:20EST qr.lisp ;; Time-stamp: <2009-12-22 22:34:28EST qr.lisp> ;; \$Id\$ (in-package :gsl) ;;; /usr/include/gsl/gsl_linalg.h ;;; FDL ;;; A general rectangular M-by-N matrix A has a ;;; QR decomposition into the product of an orthogonal ;;; M-by-M square matrix Q (where Q^T Q = I) and ;;; an M-by-N right-triangular matrix R, A = Q R. ;;; This decomposition can be used to convert the linear system A x = b ;;; into the triangular system R x = Q^T b, which can be solved by ;;; back-substitution. Another use of the QR decomposition is to ;;; compute an orthonormal basis for a set of vectors. The first N ;;; columns of Q form an orthonormal basis for the range of A, ;;; ran(A), when A has full column rank. (defmfun QR-decomposition (A &optional (tau (make-marray 'double-float :dimensions (min (dim0 A) (dim1 A))))) "gsl_linalg_QR_decomp" (((mpointer A) :pointer) ((mpointer tau) :pointer)) :inputs (A) :outputs (A tau) :documentation ; FDL "Factorize the M-by-N matrix A into the QR decomposition A = Q R. On output the diagonal and upper triangular part of the input matrix contain the matrix R. The vector tau and the columns of the lower triangular part of the matrix A contain the Householder coefficients and Householder vectors which encode the orthogonal matrix Q. The vector tau must be of length k=min(M,N). The matrix Q is related to these components by, Q = Q_k ... Q_2 Q_1 where Q_i = I - tau_i v_i v_i^T and v_i is the Householder vector v_i = (0,...,1,A(i+1,i),A(i+2,i),...,A(m,i)). This is the same storage scheme as used by lapack. The algorithm used to perform the decomposition is Householder QR (Golub & Van Loan, Matrix Computations, Algorithm 5.2.1).") (defmfun QR-solve (QR tau b &optional x-spec &aux (x (if (eq x-spec t) (make-marray 'double-float :dimensions (dimensions b)) x-spec))) ("gsl_linalg_QR_svx" "gsl_linalg_QR_solve") ((((mpointer QR) :pointer) ((mpointer tau) :pointer) ((mpointer b) :pointer)) (((mpointer QR) :pointer) ((mpointer tau) :pointer) ((mpointer b) :pointer) ((mpointer x) :pointer))) :inputs (QR tau b x) :outputs (x) :return ((or x b)) :documentation ; FDL "Solve the square system A x = b using the QR decomposition of A into (QR, tau) given by QR-decomp. The least-squares solution for rectangular systems can be found using QR-lssolve. If x-spec is NIL (default), the solution will replace b. If x-spec is T, then an array will be created and the solution returned in it. If x-spec is a marray, the solution will be returned in it. If x-spec is non-NIL, on output the solution is stored in x and b is not modified. The solution is returned from the function call.") (defmfun QR-solve-least-squares (QR tau b &optional (x (make-marray 'double-float :dimensions (dim1 QR))) (residual (make-marray 'double-float :dimensions (dim0 QR)))) "gsl_linalg_QR_lssolve" (((mpointer QR) :pointer) ((mpointer tau) :pointer) ((mpointer b) :pointer) ((mpointer x) :pointer) ((mpointer residual) :pointer)) :inputs (QR tau b) :outputs (x residual) :documentation ; FDL "The least squares solution to the overdetermined system A x = b where the matrix A has more rows than columns. The least squares solution minimizes the Euclidean norm of the residual, ||Ax - b||.The routine uses the QR decomposition of A into (QR, tau) given by #'QR-decomposition. The solution is returned in x. The residual is computed as a by-product and stored in residual.") (defmfun QR-QTvector (QR tau v) "gsl_linalg_QR_QTvec" (((mpointer QR) :pointer) ((mpointer tau) :pointer) ((mpointer v) :pointer)) :inputs (QR tau) :outputs (v) :documentation ; FDL "Apply the matrix Q^T encoded in the decomposition (QR, tau) to the vector v, storing the result Q^T v in v. The matrix multiplication is carried out directly using the encoding of the Householder vectors without needing to form the full matrix Q^T.") (defmfun QR-Qvector (QR tau v) "gsl_linalg_QR_Qvec" (((mpointer QR) :pointer) ((mpointer tau) :pointer) ((mpointer v) :pointer)) :inputs (QR tau) :outputs (v) :documentation ; FDL "Apply the matrix Q encoded in the decomposition (QR, tau) to the vector v, storing the result Q v in v. The matrix multiplication is carried out directly using the encoding of the Householder vectors without needing to form the full matrix Q.") (defmfun QR-Rsolve (QR b &optional x-spec &aux (x (if (eq x-spec t) (make-marray 'double-float :dimensions (dimensions b)) x-spec))) ("gsl_linalg_QR_Rsvx" "gsl_linalg_QR_Rsolve") ((((mpointer QR) :pointer) ((mpointer b) :pointer)) (((mpointer QR) :pointer) ((mpointer b) :pointer) ((mpointer x) :pointer))) :inputs (QR b x) :outputs (b x) :return ((or x b)) :documentation ; FDL "Solve the triangular system R x = b for x. It may be useful if the product b' = Q^T b has already been computed using QR-QTvec. If x-spec is NIL (default), the solution will replace b. If x-spec is T, then an array will be created and the solution returned in it. If x-spec is a marray, the solution will be returned in it. If x-spec is non-NIL, on output the solution is stored in x and b is not modified. The solution is returned from the function call.") (defmfun QR-unpack (QR tau &optional (Q (make-marray 'double-float :dimensions (list (dim0 QR) (dim0 QR)))) (R (make-marray 'double-float :dimensions (dimensions QR)))) "gsl_linalg_QR_unpack" (((mpointer QR) :pointer) ((mpointer tau) :pointer) ((mpointer Q) :pointer) ((mpointer R) :pointer)) :inputs (QR tau) :outputs (Q R) :documentation ; FDL "Unpack the encoded QR decomposition (QR, tau) into the matrices Q and R where Q is M-by-M and R is M-by-N.") (defmfun QR-QRsolve (Q R b &optional (x (make-marray 'double-float :dimensions (dim0 b)))) "gsl_linalg_QR_QRsolve" (((mpointer Q) :pointer) ((mpointer R) :pointer) ((mpointer b) :pointer) ((mpointer x) :pointer)) :inputs (Q R b) :outputs (x) :documentation ; FDL "Solves the system R x = Q^T b for x. It can be used when the QR decomposition of a matrix is available in unpacked form as Q, R).") (defmfun QR-update (Q R w v) "gsl_linalg_QR_update" (((mpointer Q) :pointer) ((mpointer R) :pointer) ((mpointer w) :pointer) ((mpointer v) :pointer)) :inputs (Q R w v) :outputs (w Q R) :return (Q R) :documentation ; FDL "Perform a rank-1 update w v^T of the QR decomposition (Q, R). The update is given by Q'R' = Q R + w v^T where the output matrices Q' and R' are also orthogonal and right triangular. Note that w is destroyed by the update.") (defmfun R-solve (R b &optional x-spec &aux (x (if (eq x-spec t) (make-marray 'double-float :dimensions (dimensions b)) x-spec))) ("gsl_linalg_R_svx" "gsl_linalg_R_solve") ((((mpointer R) :pointer) ((mpointer b) :pointer)) (((mpointer R) :pointer) ((mpointer b) :pointer) ((mpointer x) :pointer))) :inputs (R b x) :outputs (x) :return ((or x b)) :documentation ; FDL "Solve the triangular system R x = b in-place. On input x should contain the right-hand side b, which is replaced by the solution on output. If x-spec is NIL (default), the solution will replace b. If x-spec is T, then an array will be created and the solution returned in it. If x-spec is a marray, the solution will be returned in it. If x-spec is non-NIL, on output the solution is stored in x and b is not modified. The solution is returned from the function call.") ;;; Examples and unit test, from linalg/test.c (defun test-qr-solve-dim (matrix) "Solve the linear equation using QR with the supplied matrix and a right-hand side vector which is the reciprocal of one more than the index." (let ((dim (dim0 matrix))) (multiple-value-bind (QR tau) (QR-decomposition (copy matrix)) (QR-solve QR tau (create-rhs-vector dim) T)))) (defun test-qr-qrsolve-dim (matrix) "Solve the linear equation using QR with the supplied matrix and a right-hand side vector which is the reciprocal of one more than the index." (let ((dim (dim0 matrix))) (multiple-value-bind (QR tau) (QR-decomposition (copy matrix)) (multiple-value-bind (Q R) (QR-unpack QR tau) (QR-QRsolve Q R (create-rhs-vector dim)))))) (defun test-qr-lssolve-dim (matrix) "Solve the linear equation using QR least squares with the supplied matrix and a right-hand side vector which is the reciprocal of one more than the index. Returns the solution and the residual." (let ((dim (dim0 matrix))) (multiple-value-bind (QR tau) (QR-decomposition (copy matrix)) ;; Residual not checked. (QR-solve-least-squares QR tau (create-rhs-vector dim))))) (defun test-qr-decomp-dim (matrix) "Solve the QR decomposition with the supplied matrix and a right-hand side vector which is the reciprocal of one more than the index." (multiple-value-bind (QR tau) (QR-decomposition (copy matrix)) (multiple-value-bind (Q R) (QR-unpack QR tau) (matrix-product Q R)))) (defun test-qr-update-dim (matrix) "Test QR rank-1 update; this should return a matrix with all elements near zero." (let* ((dim0 (dim0 matrix)) (dim1 (dim1 matrix)) (u (create-matrix (lambda (i) (sin (1+ i))) dim0 nil)) (v (create-matrix (lambda (i) (+ (cos (+ 2 i)) (sin (+ 3 (expt i 2))))) dim1 nil)) (qr1 (create-matrix (lambda (i j) (+ (maref matrix i j) (* (maref u i) (maref v j)))) dim0 dim1)) (qr2 (copy matrix)) (w (make-marray 'double-float :dimensions dim0))) (multiple-value-bind (QR2 tau) (QR-decomposition qr2) (multiple-value-bind (Q2 R2) (QR-unpack QR2 tau) ;; compute w = Q^T u (matrix-product Q2 u w 1.0d0 0.0d0 :trans) (QR-update Q2 R2 w v) (matrix-product Q2 R2 qr2 1.0d0 0.0d0) (elt- qr1 qr2))))) (save-test qr ;; test_QR_solve (test-qr-solve-dim *hilb2*) (test-qr-solve-dim *hilb3*) (test-qr-solve-dim *hilb4*) (test-qr-solve-dim *hilb12*) (test-qr-solve-dim *vander2*) (test-qr-solve-dim *vander3*) (test-qr-solve-dim *vander4*) (test-qr-solve-dim *vander12*) ;; test_QR_QRsolve (test-qr-qrsolve-dim *hilb2*) (test-qr-qrsolve-dim *hilb3*) (test-qr-qrsolve-dim *hilb4*) (test-qr-qrsolve-dim *hilb12*) (test-qr-qrsolve-dim *vander2*) (test-qr-qrsolve-dim *vander3*) (test-qr-qrsolve-dim *vander4*) (test-qr-qrsolve-dim *vander12*) ;; test_QR_lssolve (test-qr-lssolve-dim *m53*) (test-qr-lssolve-dim *hilb2*) (test-qr-lssolve-dim *hilb3*) (test-qr-lssolve-dim *hilb4*) (test-qr-lssolve-dim *hilb12*) (test-qr-lssolve-dim *vander2*) (test-qr-lssolve-dim *vander3*) (test-qr-lssolve-dim *vander4*) (test-qr-lssolve-dim *vander12*) ;; test_QR_decomp (test-qr-decomp-dim *m35*) (test-qr-decomp-dim *m53*) (test-qr-decomp-dim *hilb2*) (test-qr-decomp-dim *hilb3*) (test-qr-decomp-dim *hilb4*) (test-qr-decomp-dim *hilb12*) (test-qr-decomp-dim *vander2*) (test-qr-decomp-dim *vander3*) (test-qr-decomp-dim *vander4*) (test-qr-decomp-dim *vander12*) ;; test_QR_update (test-qr-update-dim *m35*) (test-qr-update-dim *m53*) (test-qr-update-dim *hilb2*) (test-qr-update-dim *hilb3*) (test-qr-update-dim *hilb4*) (test-qr-update-dim *hilb12*) (test-qr-update-dim *vander2*) (test-qr-update-dim *vander3*) (test-qr-update-dim *vander4*) (test-qr-update-dim *vander12*))