Commons Math
  1. Commons Math
  2. MATH-1101

QR and Rank-revealing QR fail to find a least-squares solution

    Details

    • Type: Bug Bug
    • Status: Closed
    • Priority: Major Major
    • Resolution: Fixed
    • Affects Version/s: 3.2
    • Fix Version/s: 3.3
    • Labels:

      Description

      QR and RRQR (rank-revealing) algorithms fail to find a least-squares solution in some cases.

      The following code:

      final RealMatrix A = new BlockRealMatrix(3, 3);
      A.setEntry(0, 0, 1);
      A.setEntry(0, 1, 6);
      A.setEntry(0, 2, 4);
      A.setEntry(1, 0, 2);
      A.setEntry(1, 1, 4);
      A.setEntry(1, 2, -1);
      A.setEntry(2, 0, -1);
      A.setEntry(2, 1, 2);
      A.setEntry(2, 2, 5);
      final RealVector b = new ArrayRealVector(new double[]

      {5, 6, 1}

      );
      final QRDecomposition qrDecomposition = new QRDecomposition(A);
      final RRQRDecomposition rrqrDecomposition = new RRQRDecomposition(A);
      final SingularValueDecomposition svd = new SingularValueDecomposition(A);
      final RealVector xQR = qrDecomposition.getSolver().solve(b);
      System.out.printf("QR solution: %s\n", xQR.toString());
      final RealVector xRRQR = rrqrDecomposition.getSolver().solve(b);
      System.out.printf("RRSQ solution: %s\n", xRRQR.toString());
      final RealVector xSVD = svd.getSolver().solve(b);
      System.out.printf("SVD solution: %s\n", xSVD.toString());

      produces

      QR solution: {-3,575,212,378,628,897; 1,462,586,882,166,368; -1,300,077,228,592,326.5}
      RRSQ solution:

      {5,200,308,914,369,308; -2,127,399,101,332,898; 1,891,021,423,407,021}

      SVD solution:

      {0.5050344462; 1.0206677266; -0.2405935347}

      Showing that QR and RRQR algorithms fail to find the least-squares solution. This can also be verified by calculating the dot product between columns of A and A*x - b:

      // x = xQR, xRRQR or xSVD
      final RealVector r = A.operate.subtract(b);
      for (int i = 0; i < x.getDimension(); ++i)

      { final RealVector columnVector = A.getColumnVector(i); assertEquals(name, 0.0, r.dotProduct(columnVector), tolerance); }

      Only SVD method passes this test with decent tolerance (1E-14 or so).

      1. math-1101-bug.java
        2 kB
        Roman Werpachowski

        Activity

        Hide
        Roman Werpachowski added a comment -

        Attached test code for convenience.

        Show
        Roman Werpachowski added a comment - Attached test code for convenience.
        Hide
        Luc Maisonobe added a comment -

        I am not sure to understand. The matrix is exactly singular here, which is correctly identified if you pass a threshold of about 2e-16 to the QRDecomposition constructor (without passing it, the default threshold is an exact 0). With the default 0 threshold, the last diagonal element is really small (8.88e-16) and using it implies computing big values.

        As all dimensions are 3, I don't understand were you intend to have a least squares solution. What is attempted here seems to be computing a full linear solution of a singular problem.

        Show
        Luc Maisonobe added a comment - I am not sure to understand. The matrix is exactly singular here, which is correctly identified if you pass a threshold of about 2e-16 to the QRDecomposition constructor (without passing it, the default threshold is an exact 0). With the default 0 threshold, the last diagonal element is really small (8.88e-16) and using it implies computing big values. As all dimensions are 3, I don't understand were you intend to have a least squares solution. What is attempted here seems to be computing a full linear solution of a singular problem.
        Hide
        Roman Werpachowski added a comment -

        In this case the description of QRDecomposition.getSolver() is misleading, since it says:

        "Get a solver for finding the A × X = B solution in least square sense." (same for RRQRDecomposition).

        Also the documentation in http://commons.apache.org/proper/commons-math/userguide/linear.html says that the QR decomposition method can handle any matrix and return a least squares solution.

        If the QRDecomposition works better with non-zero threshold, then it may be worth considering changing the default threshold value from zero, since the current default behaviour is not what the user might expect based on the documentation.

        Finally, it is always possible to have a least squares solution of || A * x - b ||^2, because the norm is bounded from below by zero and diverges to infinity as the norm of x goes to infinity. The solution may not be unique though. Also note that the SVD solver returns a much more sensible solution, which minimizes || A * x - b ||^2.

        Show
        Roman Werpachowski added a comment - In this case the description of QRDecomposition.getSolver() is misleading, since it says: "Get a solver for finding the A × X = B solution in least square sense." (same for RRQRDecomposition). Also the documentation in http://commons.apache.org/proper/commons-math/userguide/linear.html says that the QR decomposition method can handle any matrix and return a least squares solution. If the QRDecomposition works better with non-zero threshold, then it may be worth considering changing the default threshold value from zero, since the current default behaviour is not what the user might expect based on the documentation. Finally, it is always possible to have a least squares solution of || A * x - b ||^2, because the norm is bounded from below by zero and diverges to infinity as the norm of x goes to infinity. The solution may not be unique though. Also note that the SVD solver returns a much more sensible solution, which minimizes || A * x - b ||^2.
        Hide
        Luc Maisonobe added a comment -

        Documentation has been improved as of r1570994, both in the javadoc and in the user guide.

        Concerning the default threshold to 0, the decision to keep this value was taken when discussing issue MATH-664.

        Our implementation of QR decomposition does not try to circumvent singular matrices by itself. I am not sure about this feature as it may fool the user into thinking everything is OK when in fact the matrix was singular and an approximate solution was found when the user expected an exact one.

        Show
        Luc Maisonobe added a comment - Documentation has been improved as of r1570994, both in the javadoc and in the user guide. Concerning the default threshold to 0, the decision to keep this value was taken when discussing issue MATH-664 . Our implementation of QR decomposition does not try to circumvent singular matrices by itself. I am not sure about this feature as it may fool the user into thinking everything is OK when in fact the matrix was singular and an approximate solution was found when the user expected an exact one.
        Hide
        Luc Maisonobe added a comment -

        Closing all resolved issue now available in released 3.3 version.

        Show
        Luc Maisonobe added a comment - Closing all resolved issue now available in released 3.3 version.

          People

          • Assignee:
            Unassigned
            Reporter:
            Roman Werpachowski
          • Votes:
            0 Vote for this issue
            Watchers:
            2 Start watching this issue

            Dates

            • Created:
              Updated:
              Resolved:

              Development