Mailing-List: contact issues-help@commons.apache.org; run by ezmlm
Precedence: bulk
Reply-To: issues@commons.apache.org
Date: Thu, 20 Feb 2014 13:42:20 +0000 (UTC)
From: "Luc Maisonobe (JIRA)" <jira@apache.org>
To: issues@commons.apache.org
Message-ID: <JIRA.12677221.1383393346290.72689.1392903740083@arcas>
In-Reply-To: <JIRA.12677221.1383393346290@arcas>
References: <JIRA.12677221.1383393346290@arcas>
Subject: [jira] [Comment Edited] (MATH-1053) QRDecomposition.getSolver()
 should be able to find pseudo-inverse of non-square matrices
MIME-Version: 1.0
Content-Type: text/plain; charset=utf-8
Content-Transfer-Encoding: 7bit


    [ https://issues.apache.org/jira/browse/MATH-1053?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=13906952#comment-13906952 ] 

Luc Maisonobe edited comment on MATH-1053 at 2/20/14 1:41 PM:
--------------------------------------------------------------

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
{noformat}
            [ C1 Q R1  |  C1 Q R2 ]
  C.A = [ ------------+-------------]
            [ C2 Q R1  |  C2 Q R2 ]
{noformat}
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:
{noformat}
     [ 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]
{noformat}

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?


was (Author: luc):
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
{noformat}
            [ C1 Q R1  |  C1 Q R2 ]
  C.A = [ ------------+-------------]
            [ C2 Q R1  |  C2 Q R2 ]
{noformat}
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:
{noformat}
      [ 12  -51   4    1 ]
 A = [  6  167  -68   2 ]
    [  -4   24   -41   3 ]
{noformat}

> QRDecomposition.getSolver() should be able to find pseudo-inverse of non-square matrices
> ----------------------------------------------------------------------------------------
>
>                 Key: MATH-1053
>                 URL: https://issues.apache.org/jira/browse/MATH-1053
>             Project: Commons Math
>          Issue Type: Improvement
>    Affects Versions: 3.2
>            Reporter: Sean Owen
>            Priority: Minor
>         Attachments: MATH-1053.patch
>
>
> 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.


--
This message was sent by Atlassian JIRA
(v6.1.5#6160)