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