Re: [eigen] Status of the eigen solver

On Wed, Oct 1, 2008 at 6:17 PM, Benoît Jacob <jacob@xxxxxxxxxxxxxxx> wrote:

Thanks Gael for having done this. The pseudo- prefix is fine by me.

ah, in fact they also correspond to the "real Schur decomposition":

"A = VSV^T, where A, V, S and V^T are matrices that contain real numbers only. In this case, V is an orthogonal matrix, V^T is the transpose of V, and S is a block upper triangular matrix called the real Schur form. The blocks on the diagonal of S are of size 1×1 (in which case they represent real eigenvalues) or 2×2 (in which case they are derived from complex conjugate eigenvalue pairs)."

The difference is that in our case the strictly upper part of the matrix S has been applied to the matrix V which, of course, is not unitary anymore. This transformation is applied at the end of the algorithm, so we can still access to the "true" real Schur form if needed.

So now I guess we have to convert the real form to the complex form before extracting the "true" eigenvectors. The second step is well documented everywhere, but for the conversion I cannot find anything else than the manual page of the corresponding matlab function which, of course, is useless for us !

my 2 cents,
gael.

BTW, In the latest commit-digest, some screensavers get ported to eigen2:
http://commit-digest.org/issues/2008-09-21/

Cheers,
Benoit

On Wednesday 01 October 2008 12:32:00 Gael Guennebaud wrote:
> FYI:
>
> [12:18] <CIA-20> ggael * r866561 eigen2/trunk/kdesupport/eigen2/
> (Eigen/src/QR/EigenSolver.h test/eigensolver.cpp):
> [12:18] <CIA-20> Fixes in Eigensolver:
> [12:18] <CIA-20> * eigenvectors => pseudoEigenvectors
> [12:18] <CIA-20> * added pseudoEigenvalueMatrix
> [12:18] <CIA-20> * clarify the documentation
> [12:18] <CIA-20> * added respective unit test
> [12:18] <CIA-20> Still missing: a proper eigenvectors() function.
>
> would be nice to get a proper eigenvectors() function before 2.0, anybody ?
>
> also the current implementation actually comes from: "Algol procedure hqr2,
> by Martin and Wilkinson, Handbook for Auto. Comp., Vol.ii-Linear Algebra,
> and the corresponding Fortran subroutine in EISPACK", so I guess we can
> find further information about this issue there.
>
> Cheers,
> gael.
>
> On Thu, Sep 25, 2008 at 3:15 AM, Benoît Jacob <jacob@xxxxxxxxxxxxxxx> wrote:
> > OK I see...
> >
> > This seems like a very clever idea, but I am uncomfortable with that
> > stuff being called eigenvectors and eigenvalues, since that's not exactly
> > what it is. I'm all for it as long as it has a slightly different name.
> >
> > Unless, someone explains me how this variant is what's actually useful in
> > most
> > applications.
> >
> > Personnally I don't see any urgent need for that (though it interests me
> > highly) so Gael if you're too busy you can ifdef that out for the 2.0
> > release, no problem. Somebody will eventually step up to help ;)
> >
> > I just think this must be sorted out cleanly before it enters our stable
> > API.
> >
> > Cheers,
> > Benoit
> >
> > PS we're mentioned here,
> >
> > http://www.gamedev.net/community/forums/topic.asp?topic_id=509161&whichpa
> >ge=1&#3318752
> >
> > On Thursday 25 September 2008 02:48:44 Gael Guennebaud wrote:
> > > yes, in fact this piece of code has been borrowed from Jama, and they
> > > use
> >
> > a
> >
> > > trick to keep the "eigenvectors" real.
> > >
> > > so this is their notes:
> > >
> > > If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is
> > > diagonal and the eigenvector matrix V is orthogonal. That is,
> > > the diagonal values of D are the eigenvalues, and
> > > V*V' = I, where I is the identity matrix. The columns of V
> > > represent the eigenvectors in the sense that A*V = V*D.
> > >
> > > If A is not symmetric, then the eigenvalue matrix D is block
> > > diagonal with the real eigenvalues in 1-by-1 blocks and any complex
> >
> > eigenvalues,
> >
> > > a + i*b, in 2-by-2 blocks, [a, b; -b, a]. That is, if the complex
> > > eigenvalues look like
> > > <pre>
> > >
> > > u + iv . . . . .
> > > . u - iv . . . .
> > > . . a + ib . . .
> > > . . . a - ib . .
> > > . . . . x .
> > > . . . . . y
> > > </pre>
> > > then D looks like
> > > <pre>
> > >
> > > u v . . . .
> > > -v u . . . .
> > > . . a b . .
> > > . . -b a . .
> > > . . . . x .
> > > . . . . . y
> > >
> > > This keeps V a real matrix in both symmetric and non-symmetric
> > > cases, and A*V = V*D.
> > >
> > >
> > > I guess it should be fairly easy to extract the "true" eigenvectors
> > > from the current ones, but unfortunately, I won't have time to look
> > > into it soon.
> > >
> > > gael.
> > >
> > > On Thu, Sep 25, 2008 at 2:33 AM, Timothy Hunter
> >
> > <tjhunter@xxxxxxxxxxxx>wrote:
> > > > For your information, the eigen values returned seem correct (as
> > > > checked with random matrices against matlab's computations). Also,
> > > > the computed eigen vectors have no imaginary part and are different.
> > > >
> > > > On Wed, Sep 24, 2008 at 1:47 PM, Gael Guennebaud
> > > >
> > > > <gael.guennebaud@xxxxxxxxx> wrote:
> > > > > ah, thanks for the fix. Actually it is supposed to work for non
> > > > > symmetric matrices. At least it worked a while ago. I'll check that
> > > > > when I'll have time....
> > > > >
> > > > > gael.
> > > > >
> > > > > On Wed, Sep 24, 2008 at 10:41 PM, Benoît Jacob <
> >
> > jacob@xxxxxxxxxxxxxxx>
> >
> > > > > wrote:
> > > > >> On Wednesday 24 September 2008 22:14:12 Timothy Hunter wrote:
> > > > >> > Hello,
> > > > >> > I would like to compute eigen values for non-symmetric matrices.
> > > > >> > I
> > > >
> > > > tried
> > > >
> > > > >> > to
> > > > >> > use to the EigenSolver but it does not compile because it cannot
> > > > >> > find the
> > > > >> > function .block(int, int) (see attached log).
> > > > >>
> > > > >> OK, this function was recently renamed to segment(int,int).
> > > > >> Fixed in r864468.
> > > > >>
> > > > >> > Also, the
> > > > >> > tests for non self-adjoint matrices are disabled in
> > > >
> > > > test/eigensolver.cpp
> > > >
> > > > >> > .
> > > > >> > Does it mean this class is not ready yet?
> > > > >>
> > > > >> I tried re-enabling this test, and fixed it so it compiled, and
> > > > >> the
> > > >
> > > > result
> > > >
> > > > >> is
> > > > >> that with a real, non-selfadjoint matrix, the unit-test failed at
> >
> > line
> >
> > > > >> 117.
> > > > >> So please consider it as not yet ready for now.
> > > > >>
> > > > >> Gael, if it's not yet ready maybe #ifdef it out for now.
> > > > >>
> > > > >> Cheers,
> > > > >> Benoit