Menu
▾
▴
matlisp-users
|
From: <ma...@ll...> - 2002-11-21 14:26:54
|
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))))
|
|
From: Michael A. K. <ma...@ll...> - 2002-11-22 13:36:39
|
Tunc,
> When M<N, should the matrix have rank=M or does it
> not matter?
1. Good question. The answer is, PINV as written will fail or
produce a numerically unstable result due to taking the reciprocal
of a very small singular value.
2. The solution is straight forward, let all singular values < MACH-TOL
equal zero. ==> rank(A) = r <= min(M,N). And use just those
r elements to form PINV.
3. I have an all day meeting today, but will make these changes in the
next couple of days. ....BTW what is the name of the "*MACH" in
the matlisp core....
more later
mike
|
|
From: Nicolas N. <Nic...@iw...> - 2002-11-22 14:37:23
|
"Michael A. Koerber" <ma...@ll...> writes: > Tunc, > > > When M<N, should the matrix have rank=M or does it > > not matter? Tunc, I don't know about the other members of this mailing list, but what concerns me: I would like to see discussions about new add-ons on the matlisp-user list. The traffic is not large here and it would be nice for users to learn how matlisp is extended. Nicolas. |
|
From: Raymond T. <to...@rt...> - 2002-11-22 15:10:18
|
>>>>> "Nicolas" == Nicolas Neuss <Nic...@iw...> writes:
Nicolas> I don't know about the other members of this mailing list, but what
Nicolas> concerns me: I would like to see discussions about new add-ons on the
Nicolas> matlisp-user list. The traffic is not large here and it would be nice
Nicolas> for users to learn how matlisp is extended.
If you want to know how to extend matlisp, please ask and ye shall
receive.
What do you want to discuss about add-ons? How to do them? Whether
they should be added? Useful additions?
Ray
|
|
From: Michael A. K. <ma...@ll...> - 2002-11-22 21:58:11
|
Tunc, > Can you say something about the difference between > your approach (using SVD) and the LAPACK routines called > DGELS and DGELSD. Numerical stability? DGELS will handle only full rank A, DGELSD will solve multiple right hand sides of Ax = [b1 | b2 | ... ] (each bi a column vector). An entity called PINV is not returned. If PINV is not required, then DGELSD might be a better choice. > Your question on MACH pinpoints a fundamental issue in Matlisp. > Almost all (if not all) Matlisp routines use the underlying > fortran implementation to do the work. Unless anything has > changed, the MACH is the one that fortran would use by default. > If this is not sufficient, or you need explicit access to this > number, we need to do some extra work (most likely in the > config file for matlisp). I haven't studied this enough yet, but it appears that D1MACH.LISP is used via F2CL to provide MACH information to the FORTRAN routines in, e.g., TOMS. D1MACH draws its information from the CL implimentation limits. I don't know that this is the same as that produced by underlying FORTRAN compiler's treatement of D1MACH.F. Unless it is the same, the numerically accuracy of the routines that rely on it could be impacted. This is just a concern on my part, not something that I have even looked at yet. Perhaps you or Ray can comment. mike |
|
From: Raymond T. <to...@rt...> - 2002-11-22 22:53:50
|
>>>>> "Michael" == Michael A Koerber <ma...@ll...> writes:
Michael> I haven't studied this enough yet, but it appears that D1MACH.LISP is
Michael> used via F2CL to provide MACH information to the FORTRAN routines in,
Michael> e.g., TOMS. D1MACH draws its information from the CL implimentation
Michael> limits. I don't know that this is the same as that produced by
Michael> underlying FORTRAN compiler's treatement of D1MACH.F. Unless it
Michael> is the same, the numerically accuracy of the routines that rely
Michael> on it could be impacted. This is just a concern on my part, not
Yes, this is true in general, but if LAPACK needs stuff from d1mach,
it's going to get it from the Fortran d1mach, not the Lisp d1mach.
There is a problem, however, because if there is such a thing for
LAPACK, it probably hasn't been tweaked at all so whatever those the
defaults are are getting used.
Ray
|
|
From: mak%<koe...@ll...> - 2002-11-24 13:50:29
|
Tunc,
1. Over the weekend I spent some more time on the PINV question. In
response to your question about DGELSD, the more I thought about it
the more I liked it. Since the routine accepts multiple RHS, I
suppose it can be provided an identity matrix. The result should then
be the PINV matrix. If I write something like this
(DEFUN PINV (A &OPTIONAL (B NIL) (...))
...
(LET ((B (IF (NULL B) IDENTITY-FOR-A B)))
...
(CALLING-GELSD ...) ...))
then the user can say (PINV A) to get the psuedo-inverse,
or they can say (PINV A B) to solve the multiple RHS problem.
2. The use of ?GELSD should aviod (my) concerns over differences in
D1MACH values in a given CL implimentation and the associated FORTRAN.
(BTW I found that the machine values are controlled by DLAMCH.F for
LAPACK, not D1MACH.F as I thought.)
3. Due to current time constraints, I intend to:
a. change the calling convention of PINV I posted last week to
provide this capability. The computations remain essentially
unchanged.
b. over the next few weeks, develop the ?GELSD.F wrappers and
replace the internals of PINV to use these.
any thoughts or redirections concerning this?
mike
|
|
From: Tunc S. <si...@ee...> - 2002-11-26 02:14:47
|
Hi Mike. This sounds good. Will the function have to work for A,B complex (or A real, B complex and so forth). Thanks, Tunc mak%koe...@ll... wrote: > > Tunc, > > 1. Over the weekend I spent some more time on the PINV question. In > response to your question about DGELSD, the more I thought about it > the more I liked it. Since the routine accepts multiple RHS, I > suppose it can be provided an identity matrix. The result should then > be the PINV matrix. If I write something like this > > (DEFUN PINV (A &OPTIONAL (B NIL) (...)) > ... > (LET ((B (IF (NULL B) IDENTITY-FOR-A B))) > ... > (CALLING-GELSD ...) ...)) > > then the user can say (PINV A) to get the psuedo-inverse, > or they can say (PINV A B) to solve the multiple RHS problem. > > 2. The use of ?GELSD should aviod (my) concerns over differences in > D1MACH values in a given CL implimentation and the associated FORTRAN. > (BTW I found that the machine values are controlled by DLAMCH.F for > LAPACK, not D1MACH.F as I thought.) > > 3. Due to current time constraints, I intend to: > > a. change the calling convention of PINV I posted last week to > provide this capability. The computations remain essentially > unchanged. > > b. over the next few weeks, develop the ?GELSD.F wrappers and > replace the internals of PINV to use these. > > any thoughts or redirections concerning this? > > mike > > ------------------------------------------------------- > This sf.net email is sponsored by:ThinkGeek > Welcome to geek heaven. > http://thinkgeek.com/sf > _______________________________________________ > Matlisp-users mailing list > Mat...@li... > https://lists.sourceforge.net/lists/listinfo/matlisp-users |
|
From: Michael A. K. <ma...@ll...> - 2002-11-26 13:37:35
|
Tunc, > Hi Mike. This sounds good. Will the function have to > work for A,B complex (or A real, B complex and so forth). Yes...all of those cases will be covered. mike ************************************************** Dr Michael A. Koerber Micro$oft Free Zone MIT/Lincoln Laboratory |
Thanks for helping keep SourceForge clean.
X