Re: [Lapackpp-devel] Pseudoinverse, Eigenvalues and Eigenvectors

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

Am Mittwoch, 30. M=E4rz 2005 22:26 schrieb Jacob \(Jack\) Gryn:
> I guess my vision is to ultimately have something at a more abstracted
> level, like matlab, where you could write, in C, something like C=3Dinv(A=
*B).

Ok, but that's not for lapackpp. Some of the functions in Matlab, like the =
one=20
you mentioned, are probably horribly inefficient/unstable from a numerical=
=20
point of view.

> With respect to the inverse functions, I don't see why 2 lapack functions
> can't be combined to make life simpler for people when they have the choi=
ce
> not to use it.  In my particular application, it's not run very often
> (maybe once in the entire program), so I would rather make my code easier
> to read in such a case.
>
> Ultimately, I'll leave the decision up to you, if you feel it's
> inappropriate, I'll just stick the function in my own local library
> instead, along with the pseudoinverse.

I stick to my decision that there shouldn't be additional inverse functions=
 in=20
lapackpp if they aren't in the fortran lapack.

Christian

>
> Jack
>
> > -----Original Message-----
> > From: Christian Stimming [mailto:sti...@tu...]
> > Sent: March 30, 2005 3:01 PM
> > To: Jacob (Jack) Gryn
> > Cc: lap...@li...
> > Subject: Re: [Lapackpp-devel] Pseudoinverse, Eigenvalues and Eigenvecto=
rs
> >
> > Hi,
> >
> > Am Mittwoch, 30. M=E4rz 2005 19:37 schrieb Jacob (Jack) Gryn:
> > > I'll work on LaLUInverseIP(LaGenMatDouble &A, LaVectorLongInt &PIV) a=
nd
> > > LaLUInverseIP(LaGenMatDouble &A, LaVectorLongInt &PIV, LaVectorDouble
> > > &work)
> >
> > Good.
> >
> > > How about adding
> > > LaInverse(const LaGenMatDouble &A, LaGenMatDouble &invA)
> > > LaInverseIP(LaGenMatDouble &A)
> > > For those who want to do it all in one shot?
> >
> > Sigh... do you desperately need these? I would rather prefer not to off=
er
> > these as extra functions, because it's just totally against the "spirit"
> > of a
> > numerically efficient library. In 95% of the cases people are better of
> > by using the solver functions... and already the existence of the inver=
se
> > functions will distract from that fact.
> >
> > > I'll also work on
> > >
> > > void LaEigSolve(const LaGenMatDouble &A, LaVectorDouble &eigvals,
> > > LaGenMatDouble &eigvec)
> >
> > That's good. I'm looking forward to your contributions.
> >
> > Christian
> >
> > > Jack
> > >
> > > > -----Original Message-----
> > > > From: lap...@li... [mailto:lapackpp-
> >
> > devel-
> >
> > > > ad...@li...] On Behalf Of Christian Stimming
> > > > Sent: March 30, 2005 3:55 AM
> > > > To: Jacob (Jack) Gryn
> > > > Cc: lap...@li...
> > > > Subject: Re: [Lapackpp-devel] Pseudoinverse, Eigenvalues and
> >
> > Eigenvectors
> >
> > > > Hi,
> > > >
> > > > Jacob (Jack) Gryn schrieb:
> > > > > Apparently Lapack has a function named DGETRI for obtaining the
> >
> > inverse
> >
> > > > of a
> > > >
> > > > > function.  Would it be reasonable to create a function called inv=
()
> >
> > in
> >
> > > > > LaGenMatDouble that makes use of DGETRI?  (If not somewhere else?=
).
> > > >
> > > > Ok, I didn't know of DGETRI. Yes, that's reasonable to be
> > > > implemented. However, the input to DGETRI is not any general matrix
> > > > but the input
> >
> > is
> >
> > > > the result of the DGETRF operation, which in lapack++ is available
> >
> > e.g.
> >
> > > > in LUFactorizeIP in linslv.cc and laslv.h (huh, obviously I forgot =
to
> > > > add the leading "La" to the name; sorry).
> > > >
> > > > Therefore I would suggest to add a function like
> > > > LaLUInverseIP(LaGenMatDouble &A, LaVectorLongInt &PIV) to
> > > > laslv.h/linslv.cc, where the documentation clearly says that the
> > > > input variables *must* come from LUFactorizeIP. I think that's bett=
er
> > > > than adding an inv() method to LaGenMatDouble, because that one wou=
ld
> >
> > always
> >
> > > > have to compute both dgetrf and dgetri but the fact that these are
> > > > two separate steps should be visible to lapackpp's users, IMHO.
> > > >
> > > > You can decide for yourself what to do with the WORK array of DGETRI
> > > > -
> >
> > -
> >
> > > > if this function is called several times, then you can gain a
> > > > significant speedup if the application re-uses the same WORK array
> > > > instead of allocating a new one each time (really!). Probably you c=
an
> > > > add two functions for the Inverse, one with the WORK array as
> > > > argument and the other one without, where the latter internally
> > > > allocates a new WORK array and then calls the former function.
> > > >
> > > > > In addition, I would like to build a pseudoinverse function as we=
ll
> > > > > (not implemented in lapack).  The two options would be to do it
> > > > > according to
> > > >
> > > > the
> > > >
> > > > > definition pinv(A)=3Dinv(A'A)A', or to do the following (sorry for
> > > > > the
> > > >
> > > > matlab
> > > >
> > > > > style notation) [U,D,V']=3Dsvd(A).  pinv(A)=3DVpinv(D)U', where t=
o get
> > > >
> > > > pinv(D),
> > > >
> > > > > just find the inverse of each non-zero value in the diagonal
> > > > > matrix.
> > > > >
> > > > > The second way apparently is more computationally efficient, but a
> >
> > loop
> >
> > > > with
> > > >
> > > > > if statements to get the pseudoinverse of D may take longer than
> >
> > using
> >
> > > > > lapack to calculate the pseudoinverse by definition.
> > > >
> > > > I'm still quite hesistant to add functions concerning "the inverse"
> > > > to
> >
> > a
> >
> > > > numerical library (the exception being cases like above, when there
> > > > already exists a LAPACK function for this). Really, all the numerics
> > > > guy's faces turn red and they start sweating once you tell them you
> >
> > want
> >
> > > > to calculate the inverse :-))
> > > > http://www.netlib.org/lapack/lug/node26.html#1232
> > > >
> > > > But seriously, I found for myself that the definition of "the useful
> > > > pseodo-inverse" depends much more on your actual application than y=
ou
> > > > might think, which means the calculation of the pseudoinverse should
> >
> > be
> >
> > > > done in application code instead of lapackpp code.
> > > >
> > > > The SVD way of calculating it is a good example, because when you
> > > > want to calculate pinv(D) by inverting "each non-zero value", you
> > > > have to think twice what "non-zero" actually means. Does it mean
> > > > "exactly
> >
> > equal
> >
> > > > to the floating point number 0.0"? Probably not, but rather "smaller
> >
> > in
> >
> > > > absolute value than some value epsilon". If epsilon is the machine's
> > > > precision, then it's approx. eps=3D2e-16 for doubles. But for your
> > > > application it might turn out that a much better value is something
> > > > smaller err larger, e.g. 1e-8 or 1e-4. But the whole point is that
> >
> > this
> >
> > > > really really really depends on your application. Also on
> > > > http://www.netlib.org/lapack/lug/node53.html it says that there is =
no
> > > > explicit LAPACK function for the pseudoinverse but "The effective
> >
> > rank,
> >
> > > > k, of A can be determined as the number of singular values which
> >
> > exceed
> >
> > > > a suitable threshold." So I would ask you not to implement such a
> > > > function in lapackpp.
> > > >
> > > > > Secondly, if I want to obtain the eigenvalues and eigenvectors of=
 a
> > > >
> > > > general
> > > >
> > > > > matrix, is DGEEV the correct function to implement?  Should I
> > > > > create
> >
> > a
> >
> > > > > function to add to EigSlv.cc ?
> > > >
> > > > Yes, and yes. (Or maybe dgeevx?) For dgeev/dgeevx (whichever seems
> > > > easier for you) you can add a function right away. Just go ahead.
> > > >
> > > > Thanks,
> > > >
> > > > Christian
> > > >
> > > >
> > > >
> > > >
> > > > -------------------------------------------------------
> > > > SF email is sponsored by - The IT Product Guide
> > > > Read honest & candid reviews on hundreds of IT Products from real
> >
> > users.
> >
> > > > Discover which products truly live up to the hype. Start reading no=
w.
> > > > http://ads.osdn.com/?ad_id=3D6595&alloc_id=3D14396&op=3Dclick
> > > > _______________________________________________
> > > > lapackpp-devel mailing list
> > > > lap...@li...
> > > > https://lists.sourceforge.net/lists/listinfo/lapackpp-devel
> > >
> > > -------------------------------------------------------
> > > SF email is sponsored by - The IT Product Guide
> > > Read honest & candid reviews on hundreds of IT Products from real
> > > users. Discover which products truly live up to the hype. Start readi=
ng
> > > now. http://ads.osdn.com/?ad_id=3D6595&alloc_id=3D14396&op=3Dclick
> > > _______________________________________________
> > > lapackpp-devel mailing list
> > > lap...@li...
> > > https://lists.sourceforge.net/lists/listinfo/lapackpp-devel