Re: [Lapackpp-devel] Pseudoinverse, Eigenvalues and Eigenvectors
Status: Beta
Brought to you by:
cstim
Menu
▾
▴
Re: [Lapackpp-devel] Pseudoinverse, Eigenvalues and Eigenvectors
|
From: Christian S. <sti...@tu...> - 2005-03-31 08:42:46
|
Jacob (Jack) Gryn schrieb: > 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 d= o > this, but I'm not sure if this is considered efficient for lapackpp's > purposes. No, I don't have any further ideas besides the for-loop. In such an=20 exceptional case an extra for loop in lapackpp will be fine. I don't understand why DGEEV doesn't write the eigenvalues into a=20 complex vector in the first place, but maybe that's because they don't=20 want to mix and match different data types. Whatever. In Lapackpp, both=20 real and complex matrices are available, so we return the eigenvalues in=20 complex-valued. Note, however, that this wouldn't work if someone were=20 to use lapackpp without the LA_COMPLEX_SUPPORT enabled. But that's okay,=20 IMHO. Christian PS: Sorry for not answering to your message=20 http://sourceforge.net/mailarchive/message.php?msg_id=3D11268095 -- I jus= t=20 discovered that message waiting in my lapackpp-devel inbox unanswered.=20 Fortunately you asked the question again and it's resolved now -- I hope=20 I don't overlook messages again. >=20 >=20 >>-----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 Eigenvecto= rs >> >>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. >> >>Right, that's what happens with general (non-symmetric) matrices. >> >> >>>Maybe it would be better to implement it using dgesdd? This may need = a >>>transpose of the V matrix though, which is probably bad. >> >>No, SVD is something different. >> >> >>>Alternatively, I could have it return a LaGenMatComplex for the >> >>eigenvector >> >>>component. >>> >>>Any thoughts on this? >> >>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. >> >>Christian >> >> >>>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 >=20 >=20 >=20 >=20 >=20 |
