Re: [Matlisp-users] PINV, initial release

SourceForge Headquarters 225 Broadway Suite 1600 San Diego, CA 92101 +1 (858) 422-6466

Hi Mike;  I can enter the changes.  But if you're going
to update PINV in few weeks I'd rather do it then, to avoid
doing it twice.  --Thanks, Tunc

ma...@ll... wrote:
> 
> All,
> 
> 1.  Below are proposed changes to add pseudo-inverse, PINV, to matlisp.
> 
> 2.  I've tested this on some random real and complex matrix and compared
>     the results to OCTAVE including rank deficient matricies (they match).
> 
> 3.  My intentions are to replace PINV in the next few weeks with direct
>     calls to DGELSD.F/ZGELSD.F.  I intend _no_ change to the calling or
>     return arguments.
> 
> 4.  I'd specifically like feedback on the inclusion of the
>     LEAST-SQUARES option; MTHD.  I almost removed this block of code
>     from PINV several times as being too trivial and adding, otherwise,
>     unnecessary code.  Its in the below code only because I posted a
>     similar routine last week.
> 
> 5.  Ray, Tunc: If there are no objections, or agreed to changes I need
>     to make, by the end of the week of 12/2/02, could you enter the
>     below changes in CVS?
> 
> tnx
> 
> mike
> 
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> ;;;                     CHANGES TO SYSTEM.DCL
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> diff -c ../matlisp system.dcl
> 
> *** ../matlisp/system.dcl       Tue Oct  8 19:30:39 2002
> --- system.dcl  Fri Nov 29 11:39:53 2002
> ***************
> *** 218,224 ****
>                      "mdivide"
>                      "msqrt"
>                      #-:mswindows "fft"
> !                    "geqr"))
>          (:module "special-functions"
>                 :source-pathname "matlisp:src"
>                 :binary-pathname ""
> --- 218,225 ----
>                      "mdivide"
>                      "msqrt"
>                      #-:mswindows "fft"
> !                    "geqr"
> !                    "pinv"))
>          (:module "special-functions"
>                 :source-pathname "matlisp:src"
>                 :binary-pathname ""
> 
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> ;;;                    CHANGES TO PACKAGES.LISP
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> diff -c ../matlisp packages.lisp
> 
> *** ../matlisp/packages.lisp    Tue Oct  8 19:30:39 2002
> --- packages.lisp       Fri Nov 29 11:36:07 2002
> ***************
> *** 285,290 ****
> --- 285,291 ----
>        "NUMBER-OF-ROWS"
>        "ONES"
>        "PRINT-ELEMENT"
> +      "PINV"
>        "QR"
>        "QR!"
>        "GEQR!"
> 
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> ;;;               NEW FILE FOR matlisp/src/pinv.lisp
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> (in-package "MATLISP")
> 
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> (defun pinv (a &optional (b nil) (tol nil) (mthd :svd))
>   "
>   SYNTAX
>   ======
>   (PINV A [B TOL MTHD])
> 
>   INPUT
>   -----
>   A         A Matlisp matrix, M x N
>   B         Optional.  Default NIL.  If provided, return solution to
>             A X = B.
>   TOL       Optional.  Default NIL.  Set limit on range of singular values
>             used in computation of pseudo-inverse.
>   MTHD      Optional.  Default :SVD.  Indicates the type of computation
>             to use for the pseudo-inverse.
> 
>               :LS  Force least squares interpretation
>               :SVD Force SVD based computation (default)
> 
>   OUTPUT  (VALUES PINV RANK SIGMA)
>   ------
>   PINV      Pseudo-inverse or solution to A X = B, a Matlisp matrix.
>   RANK      The effective rank of A, r = rank(A).  NIL if MTHD is :LS
>   SIGMA     The list of lenght min(M,N) of the signular values of A.
>             NIL if MTHD is :LS
> 
>   DESCRIPTION
>   ===========
>   Given the equation A X = B, solve for the value of A^+ which will
>   solve X = A^+ B.  A is assumed M x N.  X is N x L and B M x L.  The
>   solution is computed via a signular value decomposition, SVD.  Let
> 
>   A = U S V^H
> 
>   define S^+ = diag(1/s1, 1/s2, ... 1/sr, 0 ... 0) where r is rank(A)
>   Then the p-inv, A^+ = V S^+ U^H.
> 
>   The computuation of r = rank(A) is controlled by TOL.  If TOL is NIL,
>   or TOL < 0, or TOL > 1, TOL is replaced with DOUBLE-FLOAT-EPSILON.
>   All signular values <= TOL*S(1) are considered to be equal
>   to zero.  (Note: singular values have a descending order
>   S(1) >= S(2) ...S(r) > 0.)
> 
>   If A is M x M, and full rank, (M:M/ A B) can be used.
> 
>   Use MTHD :SVD when rank(A) = r < N.  For r <= M < N, this is the
>   underconstrained case for which A^+ is the minimum norm solution.
> 
>   The MTHD :LS can be used when M > N and rank(A) = N.  This is the
>   canonical least squares problem.  In this case A^+ = inv(A^H A) A^H.
>   NOTE:  NO RANK CHECKS OR SUCH ARE MADE HERE.  IF A^H A IS RANK
>   DEFICIENT THIS METHOD MAY FAIL.  USER BEWARE."
> 
>   ;; reset the TOL value as required.  This odd screening
>   ;; is for future compatibility to the ?GELSD.F routines
>   (if (or (null tol)
>           (>= tol 1)
>           (< tol 0))
>       (setf tol double-float-epsilon))
> 
>   (cond
>     ;; Force the least squares interpretation.
>    ((eq mthd :ls)
>     (let* ((ah (ctranspose a))
>            (pinv (m/ (m* ah a) ah)))
>       (values (if b (m* pinv b) pinv) nil nil)))
> 
>    ;; This is the major routine here, using the SVD method
>    ((eq mthd :svd)
>     (multiple-value-bind(up sp vptranspose info)
>         (svd a :s)                      ; use economy SVD
> 
>       (if (null info) (error "Computation of SVD failed."))
> 
>       ;; This form sets all signular values above the rank
>       ;; of A to zero
>       (let ((pinv nil)
>             (sp-values (m::store (diag sp)))
>             (rank
>              ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>              ;; This form computes the reciprocal of the singular values
>              ;; up to and including the rank, r, and returns the RANK
>              (do ((limit (* tol (matrix-ref sp 0 0)))
>                   (max-rank (min (number-of-rows sp) (number-of-cols sp)))
>                   (n 0 (incf n)))
>                  ;; terminate when all sigular values are used,
>                  ;; or fall below the limit
>                  ((or (= n max-rank)
>                       (<= (matrix-ref sp n n) limit)) n)
> 
>                (setf (matrix-ref sp n n) (/ (matrix-ref sp n n))))))
> 
>         ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>         ;; set the values of the reciprocal singular values to zero
>         ;; above the RANK of A
>         (do ((n rank (incf n))
>              (max-rank (min (number-of-rows sp) (number-of-cols sp))))
>             ((= n max-rank))
> 
>           (setf (matrix-ref sp n n) 0d0))
> 
>         ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>         ;; compute the PINV or PINV * B
>         (setf pinv (m* (ctranspose vptranspose) (m* sp (ctranspose up))))
>         (values
>          (if b (m* pinv b) pinv)
>          rank
>          sp-values))))
> 
>    (t (error "Invalid Method (~A) passed to PINV." mthd))))
> 
> -------------------------------------------------------
> This SF.net email is sponsored by: Get the new Palm Tungsten T
> handheld. Power & Color in a compact size!
> http://ads.sourceforge.net/cgi-bin/redirect.pl?palm0002en
> _______________________________________________
> Matlisp-users mailing list
> Mat...@li...
> https://lists.sourceforge.net/lists/listinfo/matlisp-users