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

[Christian S. <sti...@tu...>]

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