Commons Math
  1. Commons Math
  2. MATH-1053

QRDecomposition.getSolver() should be able to find pseudo-inverse of non-square matrices

    Details

    • Type: Improvement Improvement
    • Status: Closed
    • Priority: Minor Minor
    • Resolution: Fixed
    • Affects Version/s: 3.2
    • Fix Version/s: 3.3
    • Labels:
      None

      Description

      I don't have a complete solution to this, so don't commit this as-is, but posting in case someone can get it over the line.

      If you process a tall m x n matrix (non-square, m>n) with QRDecomposition and then call getSolver().getInverse(), you will get DimensionMismatchException. There's not a good reason the QR decomposition can't compute the least-squares solution here.

      The issue is that it tries to invert A by solving AX = I. The dimension of I has to match the row dimension of A, or m. However it's using the length of the diagonal of R, which is min(m,n), which is n when m>n.

      That patch is simple and is part of the attached patch. It also includes a test case for a tall matrix.

      However it doesn't work for a fat matrix (m<n). There's a test case for that too. It returns an n x m value but the rows for i >= m are 0 and are not computed. I'm not sure enough about the shape of the computation to be able to fix it, but it is where it's solving the triangular system Rx = y.

        Activity

        Hide
        Sean Owen added a comment -

        As a more conservative change for now, it could merely reject a tall matrix with an exception, and document this. That would probably be nicer than the current DimesnionMismatchException.

        Show
        Sean Owen added a comment - As a more conservative change for now, it could merely reject a tall matrix with an exception, and document this. That would probably be nicer than the current DimesnionMismatchException.
        Hide
        Gilles added a comment -

        This would require changing the contract of the interface.
        Not a problem, if this would indeed be better. Should be raised on the ML.

        Show
        Gilles added a comment - This would require changing the contract of the interface. Not a problem, if this would indeed be better. Should be raised on the ML.
        Hide
        Luc Maisonobe added a comment - - edited

        I have looked at the patch and it seems correct to me for tall matrices.
        Concerning fat matrices, I am not sure what you want is even possible.
        It seems to me that in the test testInvertShortWide, you computed your reference matrix by computing the decomposition of the tall transpose matrix and then transposed the result back. Then, you get a non null fourth row. The current code (with your patch applied) would compute a zero fourth row, but the three first rows would be identical.

        Regardless of this fourth row, I think the matrix can never be considered a pseudo inverse in the way you want. Here is my rationale to explain this. Consider a fat matrx A = [ A1 | A2 ], with A1 a square nxn matrix and A2 the (m-n)xn matrix containing the remaining n-m columns. QR decomposition of this matrix is A = QR with Q a nxn square matrix and R = [ R1 | R2] a fat matrix with R1 an nxn upper triangular and R2 the remaining columns.

        It is possible to find a tall matrix B such that A.B = I, but it is not possible to find a tall matrix C such that C.A = I. Here again, I will consider subscript 1 is for the square part and 2 for the remaining elements (here it will be rows).

        Block multiplication leads to A.B = Q R1 B1 + Q R2 B2 So choosing B1 by solving Q R1 B1 = I (i.e. B1 is the inverse of Q R1) and choosing B2 = 0 gives an acceptable solution. This is what the current code with your patch applied does.

        Block multiplication leads to

                    [ C1 Q R1  |  C1 Q R2 ]
          C.A = [ ------------+-------------]
                    [ C2 Q R1  |  C2 Q R2 ]
        

        Here again, we can choose C1 so the upper left part is identity (and we will in fact get C1 = B1). But then, the upper right part will be C1 Q R2 = (Q.R1)-1 Q R2 = R1-1 . R2 which only depends on A (and in fact only on the triangular parts of the decomposition as only R1 and R2 are involved), and which is in general non-zero.

        In your example, the matrices are:

             [ 12  -51   4    1 ]
         A = [  6  167  -68   2 ]
             [  -4   24   -41   3 ]
        
            [ -2450/175   -3675/175   2450/175]
        R1 = [     0      -30625/175 12250/175]
            [    0          0        6125/175]
        
          [ -150/175]
         R2 = [ -337/175]
             [ -541/175]
        

        and hence the three first elements of C.A are 23/2450, -149/6125, -541/6125, all non-zero. So C.A cannot be identity.

        So I will probably add you patch for the tall case, but generate an error for the fat case. Does this suit your needs?

        Show
        Luc Maisonobe added a comment - - edited I have looked at the patch and it seems correct to me for tall matrices. Concerning fat matrices, I am not sure what you want is even possible. It seems to me that in the test testInvertShortWide, you computed your reference matrix by computing the decomposition of the tall transpose matrix and then transposed the result back. Then, you get a non null fourth row. The current code (with your patch applied) would compute a zero fourth row, but the three first rows would be identical. Regardless of this fourth row, I think the matrix can never be considered a pseudo inverse in the way you want. Here is my rationale to explain this. Consider a fat matrx A = [ A1 | A2 ], with A1 a square nxn matrix and A2 the (m-n)xn matrix containing the remaining n-m columns. QR decomposition of this matrix is A = QR with Q a nxn square matrix and R = [ R1 | R2] a fat matrix with R1 an nxn upper triangular and R2 the remaining columns. It is possible to find a tall matrix B such that A.B = I, but it is not possible to find a tall matrix C such that C.A = I. Here again, I will consider subscript 1 is for the square part and 2 for the remaining elements (here it will be rows). Block multiplication leads to A.B = Q R1 B1 + Q R2 B2 So choosing B1 by solving Q R1 B1 = I (i.e. B1 is the inverse of Q R1) and choosing B2 = 0 gives an acceptable solution. This is what the current code with your patch applied does. Block multiplication leads to [ C1 Q R1 | C1 Q R2 ] C.A = [ ------------+-------------] [ C2 Q R1 | C2 Q R2 ] Here again, we can choose C1 so the upper left part is identity (and we will in fact get C1 = B1). But then, the upper right part will be C1 Q R2 = (Q.R1) -1 Q R2 = R1 -1 . R2 which only depends on A (and in fact only on the triangular parts of the decomposition as only R1 and R2 are involved), and which is in general non-zero. In your example, the matrices are: [ 12 -51 4 1 ] A = [ 6 167 -68 2 ] [ -4 24 -41 3 ] [ -2450/175 -3675/175 2450/175] R1 = [ 0 -30625/175 12250/175] [ 0 0 6125/175] [ -150/175] R2 = [ -337/175] [ -541/175] and hence the three first elements of C.A are 23/2450, -149/6125, -541/6125, all non-zero. So C.A cannot be identity. So I will probably add you patch for the tall case, but generate an error for the fat case. Does this suit your needs?
        Hide
        Luc Maisonobe added a comment -

        Fixed in subversion repository as of r1570566.

        Thanks for the report and for the patch.

        Show
        Luc Maisonobe added a comment - Fixed in subversion repository as of r1570566. Thanks for the report and for the patch.
        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:
            Sean Owen
          • Votes:
            0 Vote for this issue
            Watchers:
            3 Start watching this issue

            Dates

            • Created:
              Updated:
              Resolved:

              Development