# commons-issues mailing list archives

##### Site index · List index
Message view
Top
From "Luc Maisonobe (JIRA)" <j...@apache.org>
Subject [jira] [Comment Edited] (MATH-1053) QRDecomposition.getSolver() should be able to find pseudo-inverse of non-square matrices
Date Thu, 20 Feb 2014 13:42:20 GMT
```
[ 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)

```
Mime
View raw message