MPINVERT: Moore-Penrose Inverse of a Matrix

[Eddie W. Shore]

Also known as a pseudoinverse, the Moore-Penrose inverse of a matrix, denoted by A^+ (capital A with a supersubscript of a plus sign), is an inverse of matrix. Different from the "true" matrix inverse, the Moore-Penrose inverse allows for non-square matrices. Primarily, the Moore-Penrose inverses are calculated is assist in solving linear least-square equations.


Input: MPINVERT(matrix). Execute this from Home or CAS screen.

Program:
EXPORT MPINVERT(mtx)
BEGIN
// 2014-08-27 EWS
// Moore-Penrose Matrix Inverse
LOCAL r,c,d,n;
d:=SIZE(mtx);
r:=d(1);
c:=d(2);
n:=RANK(mtx);
CASE
IF n==c THEN
RETURN (TRN(mtx)*mtx)^-1*mtx;
END;
IF n==r THEN
RETURN TRN(mtx)*(mtx*TRN(mtx))^-1;
END;
DEFAULT
RETURN "No Solution Found";
END;

END;

Examples

Matrix:
[ [ 1, 2, 3 ] [ 3, 4, 0 ] ]

Moore-Penrose Inverse:
[ [ -8/229, 31/229 ] [ 6/229, 34/229 ] [ 75/229, -33/229 ] ]

Matrix:
[ [7, 4, 6, -7] [-1, 5, 3, 3] ]

Moore-Penrose Inverse: (to four decimal places)
[ [0.0489, -0.0338] [0.0194, 0.1092] [0.0360, 0.0600] [-0.0520, 0.0800] ]

[Namir]

In the case of square matrices, how close should the PM pseudo-inverse be to the actual matrix inverse?

Namir

[salvomic]

(08-29-2014 09:30 PM)Eddie W. Shore Wrote:  Also known as a pseudoinverse, the Moore-Penrose inverse of a matrix, denoted by A^+ (capital A with a supersubscript of a plus sign), is an inverse of matrix. ...

hi Eddie,
thank you!
very useful for me.
However I think that the "left" pseudo inverse should me

Code:


IF n==c THEN
RETURN (TRN(mtx)*mtx)^-1*TRN(mtx);

also you could use

Code:


r:=rowDim(mtx);
c:=colDim(mtx);
n:=RANK(mtx);

without calculate at all d:=SIZE(mtx) and so on...
Obviously you should set also
Local r,c,n;

I hope you could appreciate my little two cents

Salvo

[Han]

The SVD decomposition of a matrix can also be used to compute the pseudo-inverse of a matrix. If
\[ A = U \cdot \Sigma \cdot V^T \]
where \( \sigma_i \) are the singular values of the diagonal matrix \( \Sigma \), and \( U \) and \( V \) are unitary matrices.

The pseudo-inverse of \( A \) is
\[ A^{+} = V \cdot \Sigma^{+} \cdot U^T \]
where \( \Sigma^{+} \) is obtained by replacing the diagonal entries \( \sigma_i \) with \( 1/\sigma_i \) for \( \sigma_i > \epsilon \) and 0 otherwise and transposing. If \( A \) is non-singular, then all \( \sigma_i \) are non-zero and the pseudo-inverse is the same as the regular inverse. On the other hand, if \( A \) is singular, then only the "non-zero" (greater than \( \epsilon \)) singular values are inverted in computing \( A^{+} \).

It's a bit more powerful in that it even handles cases when \( A^TA \) or \( AA^T \) are not invertible. The SVD has lots of neat applications! What are they? Well, that's a google exercise left to the diligent reader :-)

[Han]

(09-01-2014 09:11 PM)Namir Wrote:  In the case of square matrices, how close should the PM pseudo-inverse be to the actual matrix inverse?

Namir

If A is square, and A is invertible, then they should be exactly equal. If A is invertible, then so is A^T and

\[ A^{+} = (A^TA)^{-1} A^{T} = A^{-1} (A^T)^{-1} A^T = A^{-1} \]

[salvomic]

(02-06-2015 07:38 PM)Han Wrote:  The SVD decomposition of a matrix can also be used to compute the pseudo-inverse of a matrix. If
\[ A = U \cdot \Sigma \cdot V^T \]
where \( \sigma_i \) are the singular values of the diagonal matrix \( \Sigma \), and \( U \) and \( V \) are unitary matrices.

The pseudo-inverse of \( A \) is
\[ A^{+} = V \cdot \Sigma^{+} \cdot U^T \]
where \( \Sigma^{+} \) is obtained by replacing the diagonal entries \( \sigma_i \) with \( 1/\sigma_i \) for \( \sigma_i > \epsilon \) and 0 otherwise and transposing. If \( A \) is non-singular, then all \( \sigma_i \) are non-zero and the pseudo-inverse is the same as the regular inverse. On the other hand, if \( A \) is singular, then only the "non-zero" (greater than \( \epsilon \)) singular values are inverted in computing \( A^{+} \).

It's a bit more powerful in that it even handles cases when \( A^TA \) or \( AA^T \) are not invertible. The SVD has lots of neat applications! What are they? Well, that's a google exercise left to the diligent reader :-)

Hi Han,
I'm interested in this extension, almost for the case in witch the matrix is not invertible and MPINVERT doesn't work.
I saw the SVD syntax. The guide writes "Description
Singular Value Decomposition. Factors an m n matrix into two matrices and a vector: {[[m m square orthogonal]],[[n n square orthogonal]], [real]}.".
Now I try svd([[1,2][1,2]]) and I get:
a matrix [[1/√2,0][1/√2,0]] then a vector [3.16227766..., 0] and at the end another matrix [[0.447213..., -0.894427...][0.894427..., 0.4472135955]]
What it the "real"?
Those values are U, ∑, V and ∑ is the diagonal (√10, 0). Is it right?

And therefore, how can we operate to get ∑, A+ and so on?
Eddie, then, could improve it's program for the case of singular non invertible matrix...

Another thing: in that case Prime apparently give a vector with 3 elements (matrix, vector, matrix): there is a simple way to extract one of them, i.e. to make ration approximation (with QPI)? Please, help.

Salvo

[Han]

(02-06-2015 08:40 PM)salvomic Wrote:  Hi Han,
I'm interested in this extension, almost for the case in witch the matrix is not invertible and MPINVERT doesn't work.
I saw the SVD syntax. The guide writes "Description
Singular Value Decomposition. Factors an m n matrix into two matrices and a vector: {[[m m square orthogonal]],[[n n square orthogonal]], [real]}.".
Now I try svd([[1,2][1,2]]) and I get:
a matrix [[1/√2,0][1/√2,0]] then a vector [3.16227766..., 0] and at the end another matrix [[0.447213..., -0.894427...][0.894427..., 0.4472135955]]
What it the "real"?

The vector of "real" values is a list of the singular values (i.e. the diagonal entries of \( \Sigma \).

Quote:Those values are U, ∑, V and ∑ is the diagonal (√10, 0). Is it right?

Yes. The help has been updated to correct the ordering of the output so that now it does in fact return \( \{ U, \Sigma, V \} \).

Quote:And therefore, how can we operate to get ∑, A+ and so on?
Eddie, then, could improve it's program for the case of singular non invertible matrix...

The pseudo inverse is \( V \Sigma^{-1} U^{H} \) where \( U^H \) is the Hermitian (conjugate transpose for complex matrices, or just transpose for real-valued matriced) and \( \Sigma^{-1} \) is simply the the transpose of \( \Sigma \) with the diagonal entry \( s_i \) replaced by \( \frac{1}{s_i} \). For values of \( s_i \) that are very small, (i.e. 0), then \( s_i \) is simply replaced by 0.

Quote:Another thing: in that case Prime apparently give a vector with 3 elements (matrix, vector, matrix): there is a simple way to extract one of them, i.e. to make ration approximation (with QPI)? Please, help.

You can use mat2list, and QPI accepts lists as input.

Currently, the SVD() command is not as robust. I recommend using: http://hpmuseum.org/forum/thread-4976.html

I have updated the program to also include pinv() for pseudo-inverse as well as pivoted QR factorization.

This article has a decent interpretation of the SVD (though nowhere complete, however): http://robotics.caltech.edu/~jwb/courses...pseudo.pdf

[salvomic]

thank you, Han,
I'll try soon your new program!

many thanks for the explanation of SVD syntax

Salvo