matlisp-users

[Matlisp-users] PINV, initial release

[<ma...@ll...>]

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))))

Re: [Matlisp-users] PINV, initial release

[Tunc S. <si...@ee...>]

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

Re: [Matlisp-users] PINV, initial release

[Raymond T. <to...@rt...>]

>>>>> "Tunc" == Tunc Simsek <si...@ee...> writes:

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

Sounds good to me.

A couple of comments:

o Can this be called pseudo-inverse instead of pinv?  My general rule
  has been if LAPACK has a certain useful routine, then use that name,
  otherwise pick a more useful, less abbreviated name. :-)

o Should the args to pinv be keyword args instead of optional?  If so,
  please make the mthd arg have the keyword :method instead.  Easier
  to figure out what it is. :-)

Everything else looks good.

Ray

Re: [Matlisp-users] PINV, initial release

[Jefferson P. <jp...@cs...>]

Raymond Toy wrote:

> 
> o Should the args to pinv be keyword args instead of optional?  If so,
>   please make the mthd arg have the keyword :method instead.  Easier
>   to figure out what it is. :-)

I thought that the MATLISP convention was to use optional arguments. 
This matches SVD, JOIN, etc.  Personally, I almost always prefer keyword 
arguments over optional.  Especially when there's more than one optional 
arg.

I have no real preference on the name, between PNV and pseudoinverse. 
PINV is nice because it matches matlab/octave, which makes it easy to 
find, if you're looking for it.  But pseudoinverse will be equally easy, 
I think.

J.