Re: [Toon-members] NonSymEigen for TooN, and (SVD full....)

[E. Rosten]

On Thu, 28 Jan 2010, Gabriel Nützi wrote:

> Hello guys
>
> I have adapted the SymEigen  for Non-Symetric Matrix, because there is no 
> Complex Support in TooN we have just Imaginary and real parts for the 
> EigenValues and Eigenvectors

In principle, TooN does support Vector<std::complex<double> >. I think 
some minor modification will be required, but not much.

One particular problem is that LAPACK seems to want to return the real and 
imaginary parts of the eigenvalues in different arrays. I find this a 
little odd seeing as FORTRAN supports complex types.

To make that work with TooN, we would need a new Allocator type which 
stores two pointers, and upon indexing constructs and returns a 
complex<double>. This would require that Vector also be parameterized on 
the data pointer. Perhaps the underlying allocator type could provide this 
type. Such parameterization would allow a Vector to be constructed off a 
pair<double*,double*> rather than a complex<double>*.

An intermediate solution would be to simply copy the two Vectors in to a 
single Vector<complex>. I can't see that this would be a problem since 
copying is a fast O(N) operation and eigendecomposition is a rather slow 
O(N^3) operation.

Is returning a Vector<complex<double> >  better than returning two 
vectors?

Making that change would allow one to be able to decouple the underlying 
datatype from the datatype returned by operator[]. This would allow one to 
take the conjugate transpose of a matrix cleanly.

Thoughts?


> --> NonSymEigen.h (it's not finished (I leave the documentation and add ons 
> for this class to other guys because I am not the pro here)  but would be 
> nice to have that in the later version of TooN :-). It has been tested and it 
> works.... :-) uses Lapack dgeev_ and sgeev_

We can leave out all the helper functions until someone needs them. Better 
to have missing features than to put in untested code.


> TestCode :
>
> Matrix<4,4> M;
> ...
> NonSymEigen<> myEig(M);
>
> Matrix<>a= myEig.get_evector();
> Vector<>b= myEig.get_evaluesIm();
> Vector<>c= myEig.get_evaluesRe();
>
> ====================================
>
>
> --> The SVD.h file contains an advanced version for full or not full 
> decomposition, so we dont have to Pad with zeros....

This seems to be directly supported by LAPACK, so it is probably more 
efficient this way than padding the matrix. If it is not more efficient, 
then the easiest solution is to have a separate FullSVD class which wraps 
SVD and pads the matrix.

> --> Dont know what you think about that but its handy :-)

I think it would be better to specify it using a boolean template 
parameter. This would ensure that the returned matrices are statically 
sized. This could be made compatible with existing code by renaming the 
SVD class, and making two wrapper classes, SVD and FullSVD.

That aside, I have a couple of qureies:

Even with full_compute=true, it appears that get_vertical() depends on the 
input matrix. Given that, it appeats that get_U and get_VT are missing 
some if(full_compute) code paths.

-Ed


> Thanks for the support and maintaining the class :-)
>
> Gabriel

--
(You can't go wrong with psycho-rats.)(http://mi.eng.cam.ac.uk/~er258)

/d{def}def/f{/Times s selectfont}d/s{11}d/r{roll}d f 2/m{moveto}d -1
r 230 350 m 0 1 179{ 1 index show 88 rotate 4 mul 0 rmoveto}for/s 12
    d f pop 235 420 translate 0 0 moveto 1 2 scale show showpage