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

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

Do you have any recommendations on how to convert the 2 real-values
matricies into a single complex-valued one?  I can make a for-loop to do
this, but I'm not sure if this is considered efficient for lapackpp's
purposes.

> -----Original Message-----
> From: Christian Stimming [mailto:sti...@tu...]
> Sent: March 30, 2005 3:48 PM
> To: Jacob (Jack) Gryn
> Cc: lap...@li...
> Subject: Re: [Lapackpp-devel] Pseudoinverse, Eigenvalues and =
Eigenvectors
>=20
> Am Mittwoch, 30. M=E4rz 2005 22:40 schrieb Jacob \(Jack\) Gryn:
> > With respect to void LaEigSolve, DGEEV/DGEEVX seem to return
> > eigenvalue/eigenvector pairs with the possibility that the =
eigenvalues
> > might be complex numbers.
>=20
> Right, that's what happens with general (non-symmetric) matrices.
>=20
> > Maybe it would be better to implement it using dgesdd?  This may =
need a
> > transpose of the V matrix though, which is probably bad.
>=20
> No, SVD is something different.
>=20
> > Alternatively, I could have it return a LaGenMatComplex for the
> eigenvector
> > component.
> >
> > Any thoughts on this?
>=20
> Yes, returning a LaGenMatComplex is probably better. However, the
> conversion
> from two real-valued vectors to one complex-valued will have to be =
done in
> Lapackpp.
>=20
> Christian
>=20
> >
> > Jack
> >
> > > -----Original Message-----
> > > From: Jacob (Jack) Gryn [mailto:jg...@cs...]
> > > Sent: March 30, 2005 12:37 PM
> > > To: 'Christian Stimming'
> > > Cc: 'lap...@li...'
> > > Subject: RE: [Lapackpp-devel] Pseudoinverse, Eigenvalues and
> Eigenvectors
> > >
> > > I'll work on LaLUInverseIP(LaGenMatDouble &A, LaVectorLongInt =
&PIV)
> and
> > > LaLUInverseIP(LaGenMatDouble &A, LaVectorLongInt &PIV, =
LaVectorDouble
> > > &work)
> > >
> > > How about adding
> > > LaInverse(const LaGenMatDouble &A, LaGenMatDouble &invA)
> > > LaInverseIP(LaGenMatDouble &A)
> > >
> > > For those who want to do it all in one shot?
> > >
> > > I'll also work on
> > >
> > > void LaEigSolve(const LaGenMatDouble &A, LaVectorDouble &eigvals,
> > > LaGenMatDouble &eigvec)
> > >
> > > 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