Re: [Matlisp-users] pseudoinverse

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

Jeff,

  below is your PINV code with some minor adjustments.  The need for
  the extra CTRANSPOSE was probably due to MATLISP returning V^H, not
  V as Octave does.  Also note that CTRANSPOSE of S is not needed
  since the singular values are real.

Ray,

  I added some comments to the PINV code, written by Jeff, and I'd like to 
  suggest that it be added to MATLISP.  

tnx

mike
**************************************************
Dr Michael A. Koerber    	Micro$oft Free Zone
MIT/Lincoln Laboratory	 	



(defun pinv (a &optional (flag :svd))
  "
  SYNTAX
  ======
  (PINV2 A :FLAG)

  INPUT
  -----
  A         A Matlisp matrix, M x N
  FLAG      A key indicating the type of computation to use for
            the pseudo-inverse.  Default :SVD

            :LS  Force least squares interpretation
            :SVD Force SVD based computation (default)

  OUTPUT
  ------
  PINV      Pseudo-inverse, a Matlisp matrix.

  DESCRIPTION
  ===========
  Given the equation Ax = b, solve for the value of A^+ which will
  solve x = A^+ b.  A is assumed M x N.  If A is M x M, and full rank,
  use (M:M/ A B).

  Use FLAG :LS when M > N and rank(A) = N.  This is the canonical least
  squares problem.  In this case A^+ = inv(A^H A) A^H.

  Use FLAG :SVD when rank(A) = r < N.  For M < N, this is the
  underconstrained case for which A^+ is the minium norm solution.
  To compute pinv, A^+,  in general, let A have an SVD

  A = USV^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."

  (cond
   ((eq flag :ls)
    ;; force the least squares interpretation.
    (let ((ah (ctranspose a)))
      (m* (m/ (m* ah a)) ah)))
   
   ((eq flag :svd)
    ;; use svd method
    (multiple-value-bind(up sp vptranspose info)
	(svd a :s)			; use economy SVD
      
      (if (null info) (error "Computation of SVD failed."))
      
      ;; form the reciprocal of the singular values
      (dotimes (n (min (number-of-rows sp) (number-of-cols sp)))
	(setf (matrix-ref sp n n) (/ (matrix-ref sp n n))))
      
      ;; compute the pinv
      (m* (ctranspose vptranspose) (m* sp (ctranspose up)))))

   (t (error "Invalid Flag (~A) passed to PINV." flag))))