Commons Math
  1. Commons Math
  2. MATH-297

Eigenvector computation incorrectly returning vectors of NaNs

    Details

    • Type: Bug Bug
    • Status: Closed
    • Priority: Major Major
    • Resolution: Fixed
    • Affects Version/s: 2.0
    • Fix Version/s: 2.1
    • Labels:
      None

      Description

      As reported by Axel Kramer on commons-dev, the following test case succeeds, but should fail:

      public void testEigenDecomposition() {
          double[][] m = { { 0.0, 1.0, -1.0 }, { 1.0, 1.0, 0.0 }, { -1.0,0.0, 1.0 } };
          RealMatrix rm = new Array2DRowRealMatrix(m);
          assertEquals(rm.toString(),
              "Array2DRowRealMatrix{{0.0,1.0,-1.0},{1.0,1.0,0.0},{-1.0,0.0,1.0}}");
          EigenDecompositionImpl ed = new EigenDecompositionImpl(rm,
              MathUtils.SAFE_MIN);
          RealVector rv0 = ed.getEigenvector(0);
          assertEquals(rv0.toString(), "{(NaN); (NaN); (NaN)}");
        }
      

      ed.getRealEigenvalues() returns the correct eigenvalues (2, 1, -1), but all three eigenvectors contain only NaNs.

        Activity

        Hide
        Bill Barker added a comment -

        This is now fixed as of R826550. Unfortunately, I lack karma to resolve this issue, so will have to hope that somebody else will resolve it for me.

        Show
        Bill Barker added a comment - This is now fixed as of R826550. Unfortunately, I lack karma to resolve this issue, so will have to hope that somebody else will resolve it for me.
        Hide
        Axel Kramer added a comment -

        If I'm expanding my testcase to the snippet below, I'm now getting an eigenvector with all "negative values" at index 1.
        I think this should be avoided.
        See also the solution computed by Ted Dunning on the mailing list:
        http://www.mail-archive.com/dev@commons.apache.org/msg12038.html

        double[][] m = {

        { 0.0, 1.0, -1.0 }

        ,

        { 1.0, 1.0, 0.0 }

        ,

        { -1.0, 0.0, 1.0 }

        };
        RealMatrix rm = new Array2DRowRealMatrix(m);
        assertEquals(rm.toString(),
        "Array2DRowRealMatrix{{0.0,1.0,-1.0},

        {1.0,1.0,0.0}

        ,{-1.0,0.0,1.0}}");
        EigenDecompositionImpl ed = new EigenDecompositionImpl(rm,
        MathUtils.SAFE_MIN);
        RealVector rv0 = ed.getEigenvector(0);
        RealVector rv1 = ed.getEigenvector(1);
        RealVector rv2 = ed.getEigenvector(2);
        assertEquals(rv0.toString(), "

        {0,58; 0,58; -0,58}

        ");
        assertEquals(rv1.toString(), "{-0; -0,71; -0,71}");
        assertEquals(rv2.toString(), "

        {0,82; -0,41; 0,41}

        ");

        Show
        Axel Kramer added a comment - If I'm expanding my testcase to the snippet below, I'm now getting an eigenvector with all "negative values" at index 1. I think this should be avoided. See also the solution computed by Ted Dunning on the mailing list: http://www.mail-archive.com/dev@commons.apache.org/msg12038.html double[][] m = { { 0.0, 1.0, -1.0 } , { 1.0, 1.0, 0.0 } , { -1.0, 0.0, 1.0 } }; RealMatrix rm = new Array2DRowRealMatrix(m); assertEquals(rm.toString(), "Array2DRowRealMatrix{{0.0,1.0,-1.0}, {1.0,1.0,0.0} ,{-1.0,0.0,1.0}}"); EigenDecompositionImpl ed = new EigenDecompositionImpl(rm, MathUtils.SAFE_MIN); RealVector rv0 = ed.getEigenvector(0); RealVector rv1 = ed.getEigenvector(1); RealVector rv2 = ed.getEigenvector(2); assertEquals(rv0.toString(), " {0,58; 0,58; -0,58} "); assertEquals(rv1.toString(), "{-0; -0,71; -0,71}"); assertEquals(rv2.toString(), " {0,82; -0,41; 0,41} ");
        Hide
        Bruce A Johnson added a comment - - edited

        There still seems to be a problem with the decomposition of some matrices. For example, the decomposition of the identity matrix

         {{1,0},{0,1}}

        yields the correct eigenvalues, but NaN for all the eigenvector elements.
        Crucially, the "isIncludedColumn" in the EigenDecompositionImplTest.java file always returns true (at least on my system) when the calculated eigenvectors have NaN elements, so is useless as a test for this problem.

        Also, I discovered this problem when getting NaN values doing an SVD of certain matrices (where each row has only one non-zero value). Since the SVD algorithm uses the EigenDecompositionImpl code, this seems to be a result of this current bug. (And ironically, I just told my students that one reason people love the SVD is that it essentially never fails).

        Show
        Bruce A Johnson added a comment - - edited There still seems to be a problem with the decomposition of some matrices. For example, the decomposition of the identity matrix {{1,0},{0,1}} yields the correct eigenvalues, but NaN for all the eigenvector elements. Crucially, the "isIncludedColumn" in the EigenDecompositionImplTest.java file always returns true (at least on my system) when the calculated eigenvectors have NaN elements, so is useless as a test for this problem. Also, I discovered this problem when getting NaN values doing an SVD of certain matrices (where each row has only one non-zero value). Since the SVD algorithm uses the EigenDecompositionImpl code, this seems to be a result of this current bug. (And ironically, I just told my students that one reason people love the SVD is that it essentially never fails).
        Hide
        Luc Maisonobe added a comment -

        I will look at this issue one issu MATH-308 has been solved.

        Show
        Luc Maisonobe added a comment - I will look at this issue one issu MATH-308 has been solved.
        Hide
        Bruce A Johnson added a comment -

        Thanks.

        The problem seems to be a divide by zero error with the variable diP1 in the following code, but this is probably as far as I'll get in debugging it:

           private void stationaryQuotientDifferenceWithShift(final double[] d, final double[] l,
                                                               final double lambda) {
                final int nM1 = d.length - 1;
                double si = -lambda;
                for (int i = 0, sixI = 0; i < nM1; ++i, sixI += 6) {
                    final double di   = d[i];
                    final double li   = l[i];
                    final double diP1 = di + si;
                    final double liP1 = li * di / diP1;
                    work[sixI]        = si;
                    work[sixI + 1]    = diP1;
                    work[sixI + 2]    = liP1;
                    si = li * liP1 * si - lambda;
                }
                work[6 * nM1 + 1] = d[nM1] + si;
                work[6 * nM1]     = si;
            }
        
        Show
        Bruce A Johnson added a comment - Thanks. The problem seems to be a divide by zero error with the variable diP1 in the following code, but this is probably as far as I'll get in debugging it: private void stationaryQuotientDifferenceWithShift(final double[] d, final double[] l, final double lambda) { final int nM1 = d.length - 1; double si = -lambda; for (int i = 0, sixI = 0; i < nM1; ++i, sixI += 6) { final double di = d[i]; final double li = l[i]; final double diP1 = di + si; final double liP1 = li * di / diP1; work[sixI] = si; work[sixI + 1] = diP1; work[sixI + 2] = liP1; si = li * liP1 * si - lambda; } work[6 * nM1 + 1] = d[nM1] + si; work[6 * nM1] = si; }
        Hide
        Bruce A Johnson added a comment -

        ..., and I note that what I presume to be the original Fortran code (in dlar1v.f of LAPACK) has several sections of code marked:

        • Runs a slower version of the above loop if a NaN is detected

        Hope this helps in resolving the issue,

        Bruce

        Show
        Bruce A Johnson added a comment - ..., and I note that what I presume to be the original Fortran code (in dlar1v.f of LAPACK) has several sections of code marked: Runs a slower version of the above loop if a NaN is detected Hope this helps in resolving the issue, Bruce
        Hide
        Jake Mannix added a comment -

        Doing the most trivial fix possible (dividing by the SAFE_MIN if diPl == 0) on the line which leads to NaN appears to correctly fix this for the particular case of decomposing diagonal matrices (at least for 2x2 and 3x3 cases I checked with a unit test).

        I'm not sure if this "fixes" the problem, because I'm not sure if I really grok the algorithm's implementation in this code.

        We could see better if I knew how often this pops up (any matrix whose tri-diagonal form is actually diagonal?)...

        Show
        Jake Mannix added a comment - Doing the most trivial fix possible (dividing by the SAFE_MIN if diPl == 0) on the line which leads to NaN appears to correctly fix this for the particular case of decomposing diagonal matrices (at least for 2x2 and 3x3 cases I checked with a unit test). I'm not sure if this "fixes" the problem, because I'm not sure if I really grok the algorithm's implementation in this code. We could see better if I knew how often this pops up (any matrix whose tri-diagonal form is actually diagonal?)...
        Hide
        Jake Mannix added a comment -

        Ok, well sadly it's easy to find an example which isn't fixed by just removing that one divide-by-zero:

         { {0, 1, 0}, {1, 0, 0}, {0, 0, 1} } 

        leads to perfectly reasonable eigenvalues (1, 1, -1), but NaN again rears its ugly head because findEigenVectors() assumes that, among other things, that the main diagonal does not start with a zero, and then divides by it.

        Not sure what the proper solution is to this, but a non-shifted LDL^t decomposition is a lot easier to understand to me than the other place where the NaN pops up, so maybe I can figure this one out on the plane ride down to ApacheCon tomorrow.

        Show
        Jake Mannix added a comment - Ok, well sadly it's easy to find an example which isn't fixed by just removing that one divide-by-zero: { {0, 1, 0}, {1, 0, 0}, {0, 0, 1} } leads to perfectly reasonable eigenvalues (1, 1, -1), but NaN again rears its ugly head because findEigenVectors() assumes that, among other things, that the main diagonal does not start with a zero, and then divides by it. Not sure what the proper solution is to this, but a non-shifted LDL^t decomposition is a lot easier to understand to me than the other place where the NaN pops up, so maybe I can figure this one out on the plane ride down to ApacheCon tomorrow.
        Hide
        Luc Maisonobe added a comment -

        Another step towards the solution has been checked in as of r885268.
        Just as Bruce noticed one month ago (thanks), there was an improvement in dstqds and dqds algorithms implemented in DLAR1V that are not in Dhillon's thesis.
        The problem is still not completely solved as for example in dimension 3 the decomposition of identity leads to 3 times the same
        vector (0, 0, 1) instead of giving (1, 0, 0), (0, 1, 0) and (0, 0, 1).

        Show
        Luc Maisonobe added a comment - Another step towards the solution has been checked in as of r885268. Just as Bruce noticed one month ago (thanks), there was an improvement in dstqds and dqds algorithms implemented in DLAR1V that are not in Dhillon's thesis. The problem is still not completely solved as for example in dimension 3 the decomposition of identity leads to 3 times the same vector (0, 0, 1) instead of giving (1, 0, 0), (0, 1, 0) and (0, 0, 1).
        Hide
        Luc Maisonobe added a comment -

        The original issue as been solved 2009-11-29 as of r885268.
        The remaining problem identified in the last comments has been moved in a separate JIRA issue: MATH-333

        Show
        Luc Maisonobe added a comment - The original issue as been solved 2009-11-29 as of r885268. The remaining problem identified in the last comments has been moved in a separate JIRA issue: MATH-333

          People

          • Assignee:
            Unassigned
            Reporter:
            Phil Steitz
          • Votes:
            0 Vote for this issue
            Watchers:
            2 Start watching this issue

            Dates

            • Created:
              Updated:
              Resolved:

              Development