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

SourceForge Headquarters 1320 Columbia Street Suite 310 San Diego, CA 92101 +1 (858) 422-6466

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

I guess this is another projected for another day :) =20

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 =
choice
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. =20

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.

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 =
Eigenvectors
>=20
> Hi,
>=20
> Am Mittwoch, 30. M=E4rz 2005 19:37 schrieb Jacob (Jack) Gryn:
> > I'll work on LaLUInverseIP(LaGenMatDouble &A, LaVectorLongInt &PIV) =
and
> > LaLUInverseIP(LaGenMatDouble &A, LaVectorLongInt &PIV, =
LaVectorDouble
> > &work)
>=20
> Good.
>=20
> > How about adding
> > LaInverse(const LaGenMatDouble &A, LaGenMatDouble &invA)
> > LaInverseIP(LaGenMatDouble &A)
> > For those who want to do it all in one shot?
>=20
> Sigh... do you desperately need these? I would rather prefer not to =
offer
> 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 inverse
> functions will distract from that fact.
>=20
> > I'll also work on
> >
> > void LaEigSolve(const LaGenMatDouble &A, LaVectorDouble &eigvals,
> > LaGenMatDouble &eigvec)
>=20
> That's good. I'm looking forward to your contributions.
>=20
> Christian
>=20
> >
> > 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 better =
than
> > > adding an inv() method to LaGenMatDouble, because that one would
> 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 =
can
> > > 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 =
well
> > > > (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 =
to 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 =
you
> > > 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 =
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
> >
> > -------------------------------------------------------
> > 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 =
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