Thread: [Lapackpp-devel] Cholesky Decomposition A = D'D
Status: Beta
Brought to you by:
cstim
Menu
▾
▴
lapackpp-devel
|
From: Jacob \(Jack\) G. <jg...@cs...> - 2005-03-15 21:50:08
|
I was wondering if anyone had implemented Cholesky Decomposition to obtain D or D' from A, in the equation A=D'D (with D' as the transpose of D)? I noticed a Cholesky factorization in linslv.cc, but I don't think it's what I'm looking for. The LAPACK functions I most likely need are either DPOTRF or DPBTRF Let me know if the existing functions can do what I want, if not, let me know if you have any tips as to how I can implement one of these functions into lapackpp. Thanks Jack |
|
From: Christian S. <sti...@tu...> - 2005-03-16 08:12:49
|
Hi, Jacob (Jack) Gryn schrieb: > I was wondering if anyone had implemented Cholesky Decomposition to=20 > obtain D or D=92 from A, in the equation A=3DD=92D (with D=92 as the tr= anspose=20 > of D)? I noticed a Cholesky factorization in linslv.cc, but I don=92t=20 > think it=92s what I=92m looking for. In linslv.cc the name only refers to the fact that the linear solver=20 internally uses such a decomposition, so that doesn't help you. > The LAPACK functions I most likely need are either >=20 > DPOTRF or DPBTRF Oh, wait, these two functions are already all that you need? "grep=20 dpotrf include/*.h" tells you that in fact these are available -- either=20 as inline functions in fmd.h or in the classes in spdfd.h or sybfd.h. However, I never tested any of these, so you should verify they really=20 do what they are supposed to. If there is any problem, I will happily=20 assist you in fixing them. Christian |
|
From: Jacob \(Jack\) G. <jg...@cs...> - 2005-03-23 22:15:12
|
Is there any documentation on the following matrices: LaSymmBandMatDouble LaSymmBandFactDouble LaSpdMatDouble LaSpdFactDouble As it stands, I don't really see anything to convert from a symmetric LaGenMatDouble to a LaSymmBandMatDouble or vice versa. =20 Is it possible to re-implement DPOTRF and/or DPBTRF to work with LaGenMatDouble's? Jack > -----Original Message----- > From: lap...@li... = [mailto:lapackpp-devel- > ad...@li...] On Behalf Of Christian Stimming > Sent: March 16, 2005 3:12 AM > To: Jacob (Jack) Gryn > Cc: lap...@li... > Subject: Re: [Lapackpp-devel] Cholesky Decomposition A =3D D'D >=20 > Hi, >=20 > Jacob (Jack) Gryn schrieb: > > I was wondering if anyone had implemented Cholesky Decomposition to > > obtain D or D' from A, in the equation A=3DD'D (with D' as the = transpose > > of D)? I noticed a Cholesky factorization in linslv.cc, but I don't > > think it's what I'm looking for. >=20 > In linslv.cc the name only refers to the fact that the linear solver > internally uses such a decomposition, so that doesn't help you. >=20 > > The LAPACK functions I most likely need are either > > > > DPOTRF or DPBTRF >=20 > Oh, wait, these two functions are already all that you need? "grep > dpotrf include/*.h" tells you that in fact these are available -- = either > as inline functions in fmd.h or in the classes in spdfd.h or sybfd.h. >=20 > However, I never tested any of these, so you should verify they really > do what they are supposed to. If there is any problem, I will happily > assist you in fixing them. >=20 > Christian >=20 >=20 >=20 > ------------------------------------------------------- > 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_ide95&alloc_id=14396&op=3Dick > _______________________________________________ > lapackpp-devel mailing list > lap...@li... > https://lists.sourceforge.net/lists/listinfo/lapackpp-devel |
|
From: Christian S. <sti...@tu...> - 2005-03-24 10:35:21
|
Jacob (Jack) Gryn schrieb: > Is there any documentation on the following matrices: > LaSymmBandMatDouble > LaSymmBandFactDouble > LaSpdMatDouble > LaSpdFactDouble Err... No. The *Fact* classes are mentioned in the "Lapack++ 1.0 User Guide", but that's about it. > As it stands, I don't really see anything to convert from a symmetric > LaGenMatDouble to a LaSymmBandMatDouble or vice versa. A general conversion doesn't make sense, does it? I guess you mean if and only if a general matrix happens to be a symmetric matrix, then you would like to be able to convert it to the LaSymmBandMatDouble? But "Band" also refers to a banded matrix, doesn't it? In that case, the internal storage structure is probably totally different anyways, so you have to make an element-wise copy. In theory some function for this could be added to Lapackpp, yes. > Is it possible to re-implement DPOTRF and/or DPBTRF to work with > LaGenMatDouble's? There is dgetrf, and LUFactorizeIP from laslv.h is using it. Is this what you were looking for? Christian |
|
From: Jacob \(Jack\) G. <jg...@cs...> - 2005-03-24 21:19:07
|
> -----Original Message----- > Err... No. The *Fact* classes are mentioned in the "Lapack++ 1.0 User > Guide", but that's about it. Strange, I did a google on LaSpdFactDouble and LaSymmBandFactDouble It turned up 0 results. > > > As it stands, I don't really see anything to convert from a symmetric > > LaGenMatDouble to a LaSymmBandMatDouble or vice versa. > > A general conversion doesn't make sense, does it? I guess you mean if > and only if a general matrix happens to be a symmetric matrix, then you > would like to be able to convert it to the LaSymmBandMatDouble? But > "Band" also refers to a banded matrix, doesn't it? In that case, the > internal storage structure is probably totally different anyways, so you > have to make an element-wise copy. In theory some function for this > could be added to Lapackpp, yes. Well, basically, I multiply a few LaGenMatDouble matrices together, and by construction the result is a Symmetric, Positive Definite matrix. Now, this matrix, I need to do a Cholesky factorization on. Just to test, I created a function in linslv.cc / laslv.h called void LaMatFactorize(const LaGenMatDouble &A,LaGenMatDouble& AF) based on the function in spdfd.h, but used it with a LaGenMatDouble as input. The result seems to be 'almost' correct. The function will only look at the upper triangle portion of the input matrix (it assumes symmetry), and the result is to be an upper-triangular matrix. The upper triangular portion of the resulting matrix is correct and does compare to the results obtained in matlab. The bottom triangle needs to be zeroed out. So, as a result, I have a working function, the problem is, if I manually zero out the bottom triangle, it may not be as optimal as lapack was intended to be. My questions now would be a) Is there a quick way of doing this zero-out? b) Should I submit this as a patch? c) If (b), is it ok in linslv.cc/laslv.h, or better elsewhere? Jack |
|
From: Christian S. <sti...@tu...> - 2005-03-26 10:03:33
|
Am Donnerstag, 24. M=E4rz 2005 22:18 schrieb Jacob (Jack) Gryn: > > -----Original Message----- > > Err... No. The *Fact* classes are mentioned in the "Lapack++ 1.0 User > > Guide", but that's about it. > > Strange, I did a google on > LaSpdFactDouble and > LaSymmBandFactDouble > > It turned up 0 results. My point is that the only additional (and old) documentation is the one on= =20 http://math.nist.gov/lapack++/ . Nothing else exists. > > > As it stands, I don't really see anything to convert from a symmetric > > > LaGenMatDouble to a LaSymmBandMatDouble or vice versa. > > > > A general conversion doesn't make sense, does it? I guess you mean if > > and only if a general matrix happens to be a symmetric matrix, then you > > would like to be able to convert it to the LaSymmBandMatDouble?=20 > > Well, basically, I multiply a few LaGenMatDouble matrices together, and by > construction the result is a Symmetric, Positive Definite matrix. > > Now, this matrix, I need to do a Cholesky factorization on. Sure, I understand. Of course Lapack[pp] has no means of knowing that this= =20 result happens to be symmetric, and it probably isn't even exactly symmetri= c=20 because of rounding errors (on the order of +-1e-15). So in any case you ha= ve=20 to construct a new LaSymmMatDouble class which you need to tell which part = of=20 your general matrix should be used. By the way, maybe the constructor which= =20 takes a (double*), i.e. a pointer to an array of doubles, is the right choi= ce=20 for this, and A(0,0) will return the pointer to the full array of the=20 original matrix: LaGenMatDouble A(10,20); // do something with A; then=20 LaSymmMatDouble B( A(0,0), A.size(0), A.size(1) ); but I'm not sure which part of A is used in B and which one is ignored (upp= er=20 vs. lower triangular matrix). You would have to look this up in the=20 documentation of the fortran lapack functions that are used. Again, why are you also talking about the banded matrices? > Just to test, I created a function in linslv.cc / laslv.h called > void LaMatFactorize(const LaGenMatDouble &A,LaGenMatDouble& AF) > based on the function in spdfd.h, but used it with a LaGenMatDouble as > input. > > The result seems to be 'almost' correct. The function will only look at > the upper triangle portion of the input matrix (it assumes symmetry), and > the result is to be an upper-triangular matrix. > > The upper triangular portion of the resulting matrix is correct and does > compare to the results obtained in matlab. The bottom triangle needs to = be > zeroed out. That results sounds 'exactly' correct to me, because that's the way how=20 fortran lapack stores symmetric matrices internally. One half of them is=20 simply disregarded, so it doesn't even need to be zeroed out. I guess you=20 refer to what you retrieve when you print it to an std::ostream? In that=20 case, the respective operator would have to treat this similar as the lapac= k=20 functions, i.e. ignoring one half and printing the full matrix solely based= =20 on the other half. > So, as a result, I have a working function, the problem is, if I manually > zero out the bottom triangle, it may not be as optimal as lapack was > intended to be. > > My questions now would be > a) Is there a quick way of doing this zero-out? Again, that part is simply ignored. It shouldn't even be zeroed out, simply= =20 ignored. > b) Should I submit this as a patch? No. See above. If you need this special case for yourself, then that's clea= ry=20 well placed in your own code. Thanks for asking, anyway. It's nice to see that the library is actually be= ing=20 used :-)) Christian |
|
From: Jacob \(Jack\) G. <jg...@cs...> - 2005-03-24 21:54:25
|
Just wondering if there's a fast way to do the following with lapack++ (if not, then in LAPACK) to do the following functions: a) Pseudoinverse (maybe I should do inv(A'A)A', but inv() needs to be done via solving the linear equation AA'=I (is there a faster way to do this?) b) Get eigenvectors and eigenvalues of a general matrix? (or even of a symmetric matrix)? Let me know Thanks Jack |
|
From: Jacob \(Jack\) G. <jg...@cs...> - 2005-03-29 22:25:36
|
Hi, 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?). 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)=inv(A'A)A', or to do the following (sorry for the matlab style notation) [U,D,V']=svd(A). pinv(A)=Vpinv(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. Any thoughts on the above? 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 ? Let me know Jack > -----Original Message----- > From: lap...@li... [mailto:lapackpp-devel- > ad...@li...] On Behalf Of Jacob (Jack) Gryn > Sent: March 24, 2005 4:54 PM > To: lap...@li... > Subject: [Lapackpp-devel] Pseudoinverse, Eigenvalues and Eigenvectors > > Just wondering if there's a fast way to do the following with lapack++ (if > not, then in LAPACK) to do the following functions: > > a) Pseudoinverse (maybe I should do inv(A'A)A', but inv() needs to be done > via solving the linear equation AA'=I (is there a faster way to do this?) > > b) Get eigenvectors and eigenvalues of a general matrix? (or even of a > symmetric matrix)? > > Let me know > > Thanks > > Jack > > > > > ------------------------------------------------------- > 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=6595&alloc_id=14396&op=click > _______________________________________________ > lapackpp-devel mailing list > lap...@li... > https://lists.sourceforge.net/lists/listinfo/lapackpp-devel |
|
From: Christian S. <sti...@tu...> - 2005-03-30 08:55:22
|
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)=inv(A'A)A', or to do the following (sorry for the matlab > style notation) [U,D,V']=svd(A). pinv(A)=Vpinv(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=2e-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 |
|
From: Jacob \(Jack\) G. <jg...@cs...> - 2005-03-30 17:37:37
|
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)=inv(A'A)A', or to do the following (sorry for the > matlab > > style notation) [U,D,V']=svd(A). pinv(A)=Vpinv(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=2e-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=6595&alloc_id=14396&op=click > _______________________________________________ > lapackpp-devel mailing list > lap...@li... > https://lists.sourceforge.net/lists/listinfo/lapackpp-devel |
|
From: Christian S. <sti...@tu...> - 2005-03-30 19:59:55
|
Hi, 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) 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 offer=20 these as extra functions, because it's just totally against the "spirit" of= a=20 numerically efficient library. In 95% of the cases people are better of by= =20 using the solver functions... and already the existence of the inverse=20 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 Eigenvecto= rs > > > > Hi, > > > > Jacob (Jack) Gryn schrieb: > > > Apparently Lapack has a function named DGETRI for obtaining the inver= se > > > > 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 lo= op > > > > 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 |
|
From: Jacob \(Jack\) G. <jg...@cs...> - 2005-03-30 20:26:38
|
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 |
×
Want the latest updates on software, tech news, and AI?
Get latest updates about software, tech news, and AI from SourceForge directly in your inbox once a month.
Thanks for helping keep SourceForge clean.
X