Commons Math
  1. Commons Math
  2. MATH-1045

EigenDecomposition.Solver should consider tiny values 0 for purposes of determining singularity

    Details

    • Type: Bug Bug
    • Status: Resolved
    • Priority: Minor Minor
    • Resolution: Fixed
    • Affects Version/s: 3.2
    • Fix Version/s: 3.3

      Description

      EigenDecomposition.Solver tests for singularity by comparing eigenvalues to 0 for exact equality. Elsewhere in the class and in the code, of course, very small values are considered 0. This causes the solver to consider some singular matrices as non-singular.

      The patch here includes a test as well showing the behavior – the matrix is clearly singular but isn't considered as such since one eigenvalue are ~1e-14 rather than exactly 0.

      (What I am not sure of is whether we should really be evaluating the norm of the imaginary eigenvalues rather than real/imag components separately. But the javadoc says the solver only supports real eigenvalues anyhow, so it's kind of moot since imag=0 for all eigenvalues.)

      1. MATH-1045.patch
        2 kB
        Sean Owen
      2. MATH-1045.patch
        3 kB
        Sean Owen
      3. MATH-1045_2.patch
        2 kB
        Sean Owen

        Issue Links

          Activity

          Hide
          Gilles added a comment -

          Isn't the code in this class (and others similarly) supposed to work for a matrix with very small entries too? I mean that, if all eigenvalues are of the order of, say, EPSILON / 10, should the matrix be considered singular right away?

          Show
          Gilles added a comment - Isn't the code in this class (and others similarly) supposed to work for a matrix with very small entries too? I mean that, if all eigenvalues are of the order of, say, EPSILON / 10, should the matrix be considered singular right away?
          Hide
          Sean Owen added a comment -

          That's a good point. If you make the example matrix non-singular, but then divide elements by 1e12, it will report it as singular. This seems wrong. On the other hand it seems a bit undesirable to return an 'inverse' in this case – it's dominated by the inverse of that tiny eigenvalue, which is huge, and the result is pretty unreliable.

          I'm a bit out of my depth here but I wonder if it's more reasonable to examine the eigenvalues in sorted order and examine ratio of one to the next. When that ratio is below epsilon it makes more sense to declare it "0".

          I could also see this being a case of "caller beware". That's the more conservative thing here.

          Show
          Sean Owen added a comment - That's a good point. If you make the example matrix non-singular, but then divide elements by 1e12, it will report it as singular. This seems wrong. On the other hand it seems a bit undesirable to return an 'inverse' in this case – it's dominated by the inverse of that tiny eigenvalue, which is huge, and the result is pretty unreliable. I'm a bit out of my depth here but I wonder if it's more reasonable to examine the eigenvalues in sorted order and examine ratio of one to the next. When that ratio is below epsilon it makes more sense to declare it "0". I could also see this being a case of "caller beware". That's the more conservative thing here.
          Hide
          Gilles added a comment -

          On the other hand it seems a bit undesirable to return an 'inverse' in this case – it's dominated by the inverse of that tiny eigenvalue, which is huge, and the result is pretty unreliable.

          This could be case-dependent and the code should perhaps be able to detect and accept input that can return a reliable result. In r1534709, I've committed an example that seems to work, even as the eigenvalues are quite small indeed.

          it's more reasonable to examine the eigenvalues in sorted order and examine ratio

          That's an interesting idea.
          Could you try and see whether it would let the new test pass, while intercepting the singular matrix of your test?

          Show
          Gilles added a comment - On the other hand it seems a bit undesirable to return an 'inverse' in this case – it's dominated by the inverse of that tiny eigenvalue, which is huge, and the result is pretty unreliable. This could be case-dependent and the code should perhaps be able to detect and accept input that can return a reliable result. In r1534709, I've committed an example that seems to work, even as the eigenvalues are quite small indeed. it's more reasonable to examine the eigenvalues in sorted order and examine ratio That's an interesting idea. Could you try and see whether it would let the new test pass, while intercepting the singular matrix of your test?
          Hide
          Sean Owen added a comment -

          On a little more research, it seems the thing to do is look at the ratio with the largest eigenvalue. Attached is a new patch where your new test and my (2) new tests all pass.

          Comments welcome from linear algebra experts but I think this is at least as principled as the existing code.

          We could also let the user specify the zero threshold as in the QRDecomposition class.

          Show
          Sean Owen added a comment - On a little more research, it seems the thing to do is look at the ratio with the largest eigenvalue. Attached is a new patch where your new test and my (2) new tests all pass. Comments welcome from linear algebra experts but I think this is at least as principled as the existing code. We could also let the user specify the zero threshold as in the QRDecomposition class.
          Hide
          Gilles added a comment -

          We could also let the user specify the zero threshold as in the QRDecomposition class.

          That would be best, I also think.
          However, there is a practical problem in that there is currently a (deprecated) constructor with the required signature.
          Could you raise this issue on the "dev" ML, and ask confirmation on how to proceed? I seem to recall that such a (functionally non-compatible) change would now be acceptable, even in a minor release.

          Show
          Gilles added a comment - We could also let the user specify the zero threshold as in the QRDecomposition class. That would be best, I also think. However, there is a practical problem in that there is currently a (deprecated) constructor with the required signature. Could you raise this issue on the "dev" ML, and ask confirmation on how to proceed? I seem to recall that such a (functionally non-compatible) change would now be acceptable, even in a minor release.
          Hide
          Thomas Neidhart added a comment -

          I just realized that in the case of a non-symmetric matrix, the eigenvalues are not sorted by value.
          The reason is that for the two cases, symmetric and non-symmetric there are completely different ways to do the decomposition.

          So you should not rely on that a particular order of the eigenvalues before we change that first.

          Show
          Thomas Neidhart added a comment - I just realized that in the case of a non-symmetric matrix, the eigenvalues are not sorted by value. The reason is that for the two cases, symmetric and non-symmetric there are completely different ways to do the decomposition. So you should not rely on that a particular order of the eigenvalues before we change that first.
          Hide
          Sean Owen added a comment -

          These are good points, as well as the comments from the thread on commons-dev – Ted in particular notes that you can use the threshold in the decomposition itself to simply stop computing eigenvalues when they get small.

          Now, these more advanced changes are near the limit of my ability and I am not sure I feel confident making them. I propose these additional changes be considered in another issue: (maybe) moving a threshold parameter, maybe modifying the decomposition.

          The patch here does I think represent a small, distinct positive change, in that it employs a reasonable test for singularity after the fact.

          Show
          Sean Owen added a comment - These are good points, as well as the comments from the thread on commons-dev – Ted in particular notes that you can use the threshold in the decomposition itself to simply stop computing eigenvalues when they get small. Now, these more advanced changes are near the limit of my ability and I am not sure I feel confident making them. I propose these additional changes be considered in another issue: (maybe) moving a threshold parameter, maybe modifying the decomposition. The patch here does I think represent a small, distinct positive change, in that it employs a reasonable test for singularity after the fact.
          Hide
          Thomas Neidhart added a comment -

          Sure np.

          I think we are just collecting all relevant information in order to decide how to proceed. The necessary changes can then be done, e.g. by myself or another maintainer.

          Show
          Thomas Neidhart added a comment - Sure np. I think we are just collecting all relevant information in order to decide how to proceed. The necessary changes can then be done, e.g. by myself or another maintainer.
          Hide
          Gilles added a comment -

          I propose these additional changes be considered in another issue

          I agree.

          Show
          Gilles added a comment - I propose these additional changes be considered in another issue I agree.
          Hide
          Gilles added a comment -

          Applied in revision 1536766.

          I created MATH-1049 for discussing further improvements.

          Show
          Gilles added a comment - Applied in revision 1536766. I created MATH-1049 for discussing further improvements.
          Hide
          Thomas Neidhart added a comment -

          We need to change the behavior for non-symmetric matrices as the eigenvalues are not sorted in this case.
          This patch relies on a descending sort order to determine if the decomposed matrix is singular, so this may fail in such a case.

          I will create a separate issue for this as it makes sense to always sort the eigenvalues imho.

          Show
          Thomas Neidhart added a comment - We need to change the behavior for non-symmetric matrices as the eigenvalues are not sorted in this case. This patch relies on a descending sort order to determine if the decomposed matrix is singular, so this may fail in such a case. I will create a separate issue for this as it makes sense to always sort the eigenvalues imho.
          Hide
          Gilles added a comment -

          [...] this may fail [...]

          For this issue, I could add a loop in order to find the one with largest absolute value. WDYT?
          But I have no idea how to construct a matrix for a unit test that would exhibit the problem.

          Show
          Gilles added a comment - [...] this may fail [...] For this issue, I could add a loop in order to find the one with largest absolute value. WDYT? But I have no idea how to construct a matrix for a unit test that would exhibit the problem.
          Hide
          Sean Owen added a comment -

          Yes, this is a good point. It's safest to find the largest eigenvalue (by absolute value) with a loop I think.

          The final matrix in testUnsymmetric(), which is unsymmetric, shows this.

          The symmetric matrix in testSquareRootNonPositiveDefinite() also shows this – the last eigenvalue is the most negative, but is the largest in absolute value.

          Show
          Sean Owen added a comment - Yes, this is a good point. It's safest to find the largest eigenvalue (by absolute value) with a loop I think. The final matrix in testUnsymmetric(), which is unsymmetric, shows this. The symmetric matrix in testSquareRootNonPositiveDefinite() also shows this – the last eigenvalue is the most negative, but is the largest in absolute value.
          Hide
          Gilles added a comment -

          Would you mind creating a patch?

          Show
          Gilles added a comment - Would you mind creating a patch?
          Hide
          Sean Owen added a comment -

          Ah right, I should have said that these show out-of-order eigenvalues but that's not what we need to check. We need one that puts a very small eigenvalue first. That's easy to generate as something like

          [ d 0 ]
          [ 1 1 ]

          for tiny d. Patch attached.

          Show
          Sean Owen added a comment - Ah right, I should have said that these show out-of-order eigenvalues but that's not what we need to check. We need one that puts a very small eigenvalue first. That's easy to generate as something like [ d 0 ] [ 1 1 ] for tiny d. Patch attached.
          Hide
          Gilles added a comment -

          Thank you. Committed in revision 1537099.

          Show
          Gilles added a comment - Thank you. Committed in revision 1537099.
          Hide
          Sean Owen added a comment -

          This is a separate issue, but so minor not sure if it merits another JIRA. While looking at this code I noticed this loop at EigenDecomposition:945 that does nothing:

          // Vectors of isolated roots
          for (int i = 0; i < n; i++) {
          if (i < 0 | i > n - 1) {
          for (int j = i; j < n; j++)

          { matrixP[i][j] = matrixT[i][j]; }

          }
          }

          The 'if' can never be true. (Not to mention non-short-circuit boolean op there.)

          Show
          Sean Owen added a comment - This is a separate issue, but so minor not sure if it merits another JIRA. While looking at this code I noticed this loop at EigenDecomposition:945 that does nothing: // Vectors of isolated roots for (int i = 0; i < n; i++) { if (i < 0 | i > n - 1) { for (int j = i; j < n; j++) { matrixP[i][j] = matrixT[i][j]; } } } The 'if' can never be true. (Not to mention non-short-circuit boolean op there.)
          Hide
          Gilles added a comment -

          issue [...] so minor not sure if it merits another JIRA

          That's certainly worth a report!
          Thanks.

          Show
          Gilles added a comment - issue [...] so minor not sure if it merits another JIRA That's certainly worth a report! Thanks.
          Hide
          Thomas Neidhart added a comment -

          This is an artifact from the original Jama source code.

          There was a similar code construct, which has been removed when the code has been translated, but this one remained. Imho it is safe to remove this part also.

          While checking this I have seen there was a new release of Jama last November with a bugfix for possible infinite loops. Need to check if our code is also affected, but most likely.

          Show
          Thomas Neidhart added a comment - This is an artifact from the original Jama source code. There was a similar code construct, which has been removed when the code has been translated, but this one remained. Imho it is safe to remove this part also. While checking this I have seen there was a new release of Jama last November with a bugfix for possible infinite loops. Need to check if our code is also affected, but most likely.
          Hide
          Thomas Neidhart added a comment -

          Removed the spurious code fragment in r1537616.

          Created and fixed also MATH-1051 to port a bugfix from Jama-1.0.3 to CM.

          Show
          Thomas Neidhart added a comment - Removed the spurious code fragment in r1537616. Created and fixed also MATH-1051 to port a bugfix from Jama-1.0.3 to CM.

            People

            • Assignee:
              Unassigned
              Reporter:
              Sean Owen
            • Votes:
              0 Vote for this issue
              Watchers:
              2 Start watching this issue

              Dates

              • Created:
                Updated:
                Resolved:

                Development