#146 failure in MatrixOps::eigens in CVS

closed-out-of-date
core (120)
6
2010-08-01
2006-06-13
No

Platform: Linux/Ubuntu/Dapper, gligc 2.3.6-0ubuntu20

The eigens function from 2.4.2 used to work on the
below symmetric marix. Now, the cvs version fails.
The below inlined perl code illustrates the error.
Other 8x8 matrixes fail silently. The
PDL::Slatec::eigsys function seems okay, so that is the
workaround. This is the error message it prints before
it prints eigenvalues and -vectors full of NaNs.

"Failure in hqr2 function. Do not trust the given
eigenvectors and -values"

I have attached a patch for t/matrixops.t which will
catch the error.

#!/usr/bin/perl -w

use strict;
use PDL;
use PDL::MatrixOps qw(eigens eigens_sym);

my @matrix;
while (<DATA>) {
my @row = split " ";
push @matrix, [@row];
}
my $m = pdl(@matrix);

print $m;
my $mt = $m->transpose();
print $mt - $m;
my ($e,$v) = eigens_sym($m);
print "eigenvalues: $e\n";
print "eigenvectors: $v\n";
($e,$v) = eigens($m);
print "eigenvalues: $e\n";
print "eigenvectors: $v\n";

__DATA__
1.638 2.153 1.482 1.695 -0.557 -2.443 -0.71 1.983
2.153 3.596 2.461 2.436 -0.591 -3.711 -0.493 2.434
1.482 2.461 2.5 2.834 -0.665 -2.621 0.248 1.738
1.695 2.436 2.834 4.704 -0.629 -2.913 0.576 2.471
-0.557 -0.591 -0.665 -0.629 19 0.896 8.622 -0.254
-2.443 -3.711 -2.621 -2.913 0.896 5.856 1.357 -2.915
-0.71 -0.493 0.248 0.576 8.622 1.357 20.8 -0.622
1.983 2.434 1.738 2.471 -0.254 -2.915 -0.622 3.214

Discussion

  • Bill Coffman

    Bill Coffman - 2006-06-21

    patch to test that catches eigens error (correction)

     
  • Bill Coffman

    Bill Coffman - 2006-06-29

    Logged In: YES
    user_id=1538801

    Further research shows that a library called ssl for "Small
    Scientific Library" was added to the cvs tree after the
    2.4.2 release.

    That subdirectory has been removed now, but the program it
    used for the eigens function remains. My tests indicate
    that this function does not converge to a solution for many
    examples that normally do. The 'eigens_sym' function, for
    symmetric matrices only, does what the old eigens function
    used to do. The new eigens is supposed to work on
    asymmetric matrices that have real eigenvalues.
    Theoretically a gain, but my analysis shows that the C code
    itself is the problem. Also, I found some errors in the
    matrixops.pd code too. Does this work: int f(int n) {
    int x[n]; x[n-1]=foo; ... Answer: *NO*.

    A workaround is to take note if the input matrix is
    symmetric, and call the older function if so. I suppose I
    can put together a fix if anyone is interested in testing
    and applying the patch. This would take the eigens function
    back to the capabilities of the 2.4.2 release.

     
  • Bill Coffman

    Bill Coffman - 2006-07-03
    • priority: 5 --> 6
    • assigned_to: nobody --> zowie
    • status: open --> pending
     
  • Bill Coffman

    Bill Coffman - 2006-07-03

    Logged In: YES
    user_id=1538801

    Test has been checked in. Code has not, but is attached, as
    well as sent to pdl-porters email. I have attached the code
    as a patch. Also, if given the word, I will check in, as I
    now have cvs write privs.

     
  • Bill Coffman

    Bill Coffman - 2006-07-03

    fix

     
  • Bill Coffman

    Bill Coffman - 2006-07-17

    Logged In: YES
    user_id=1538801

    Changing to open, since it was pending but pending what?
    Decisions have to be made about support for eigens.

    The latest changes made the functionality move back to that
    of the 2.4.2 version, but with a barf if someone attempted
    to take the eigens of an asymmetric matrix. This is due to
    deep problems with the "small scientific library" version of
    eigens, which seems to be hopelessly broken.

    To fix this bug correctly, a decision needs to be made about
    what to do when complex eigenvlaues are returned. Also, the
    simpler case of an asymmetric matrix that does not return
    complex values.

    I propose, a new function called eigens_complex that accepts
    a complex matrix and returns complex values. Also,
    eigens_real that accepts a real matrix and fails if any of
    the eigenvalues are complex. Finally an eigens_real_complex
    which takes a real matrix and returns complex numbers for
    eigenvalue/vectors. eigens should default to eigens_real.

    Okay, feel free to change the function names.

    Oh yeah, eigens_sym could take a vector with the lower
    triangular matrix in it. That's what it processes anyway,
    and many people will start with vectors anyway. Maybe
    introduce eigens_sym_lt and eigens_sym_ut.

     
  • Bill Coffman

    Bill Coffman - 2006-07-17
    • status: pending --> open
     
  • Chris Marshall

    Chris Marshall - 2006-07-18

    Logged In: YES
    user_id=44920

    I suggest fixing the cvs eigens code (and documentation)
    to work with symmetric matrices in an upcoming point
    release (basically revert to 2.4.2 documented behavior).

    Let's defer full/correct/improved eigen support for the
    next major release. I would like to see the current cvs
    with major outstanding bugs fixed as a 2.4.3 release by
    the middle of August.

    More drastic refactoring and additions would push out the
    interim release date. I think the existing cvs is
    significantly improved over 2.4.2 and worth more general
    release. When the answer to user queries is "try the CVS
    version" it is probably time to get a release out the
    door, even if only a "bugfix release".

     
  • Craig DeForest

    Craig DeForest - 2008-07-02

    Logged In: YES
    user_id=20200
    Originator: NO

    Yikes! I finally did some more systematic testing with the SSL eigens algorithm. While spot-checks with constructed matrices seem to work OK where it yields an answer, the algorithm does not coverge frequently at all. In the interest of speed (of release of 2.4.3) I'm going to leave it in switcheroo mode (using eigens_sym where possible, and the SSL where not) except that, if the GSL is available, I'll have it call GSL. That is an awful hack, but it should go away if we actually get off our duffs and integrate the GSL more thoroughly for PDL 2.5.

     
  • Craig DeForest

    Craig DeForest - 2008-07-07

    Logged In: YES
    user_id=20200
    Originator: NO

    I set up the simple switch advocated back in 2006, with a dire one-time warning message talking about the sketchy algorithm and promising to fix it for 2.5. It seems that a great deal will be overhauled for 2.5 if enough of us can get off our duffs. Leaving this open as it still needs fixing, but it's resolved "enough" for 2.4.4.

     
  • Chris Marshall

    Chris Marshall - 2009-04-29

    I was reviewing and consolidating various feature requests for better linear algebra support in PDL and decided the easiest way forward would be to extend our support for the GSL. Taking another look at the GSL reference, that would also address this problem with PDL. The consolidated GSL features request is at:

    http://sourceforge.net/tracker/?func=detail&aid=2784029&group_id=612&atid=350612

    I recommend we skip piecemeal work on various numeric issues, go with GSL as a basic capability---like FFTW, and move on.

     
  • Chris Marshall

    Chris Marshall - 2010-08-01
    • status: open --> closed-out-of-date
     
  • Chris Marshall

    Chris Marshall - 2010-08-01

    This bug is no longer valid and was specifically
    resolved for the PDL-2.4.4 and higher releases.
    We still have feature request tickets open for
    improved matrix computation support in PDL
    via GSL or other alternative library.

     

Log in to post a comment.

Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:

JavaScript is required for this form.





No, thanks