Moore–Penrose Inverse Simplified!

(Estimation of pseudo-inverse of non-square matrices (Ax=b))

Mukesh Kumar
Accredian
6 min readJun 21, 2022

--

Preface

In linear algebra, the Moore–Penrose inverse is the most widely known generalization of the inverse matrix. The Moore–Penrose inverse is sometimes also known by the name of pseudo-inverse or generalized inverse. E. H. Moore in 1920, Arne Bjerhammar in 1951, and Roger Penrose in 1955 described it independently in their research, and researchers and developers utilize it widely in solving real-life problems.

It helps estimate the best fit solution (least squares) to a system of linear equations (Ax=b).

In this article, I’ll walk you through the working of pseudo-inverse manually along with its python implementation. So, let’s dive in…

Derivation of Pseudo-inverse

Note: Skip the derivation and use the final result only if getting maths heavy.

To estimate the solution of pseudo-inverse, Singular Value Decomposition (economy variant) comes into the picture. We know that Economy SVD is

We know that the system of equation is

Replace A with the economy variant SVD

Multiply both sides by the following equation

Our equation will become:

We can re-write the above equation as

which is called Moore-Penrose Inverse or Pseudoinverse.

Manual Example of Estimating Pseudo-inverse

Let’s say we have the following matrix of 4 X 2 in size and we want to compute its pseudo-inverse.

If you are interested in how to compute the SVD of any matrix, do not worry. For now, try using the online calculator to decompose the matrix. In the future, I will also explain its decomposition. Our current objective is to find the pseudo-inverse because that’s what we want. Our decomposition matrices are:

What’s left? just using the derived formula and getting the pseudo-inverse of the matrix A. Now transposing a matrix is simple, but taking the inverse is not. So, if you don’t know how to take the inverse of a matrix, not a problem. Use the online calculator to find the inverse of a matrix. Once we find the inverse of the matrix and transpose the U and V matrices, we can proceed with estimating the pseudo-inverse.

Using the Moore-Penrose inverse equation:

We get pseudo-inverse, which is exactly what we want. Easy, right?

Python Implementation of Pseudo-inverse

Next, we will see how to implement pseudo-inverse in python.

Input:# Import Numpy library
import numpy as np

We will create a NumPy array (2d) which will contain some values in it.

Input:arr = np.array([[2, 4], [1, 3], [0, 0], [0, 0]])
arr
Output:array([[2, 4],
[1, 3],
[0, 0],
[0, 0]])

We will decompose our input matrix into U, Sigma, and V transpose as follows:

Input:U, S, VT = np.linalg.svd(a=arr)print(U)print(S)print(VT)Output:[[-0.81741556 -0.57604844  0.          0.        ]  
[-0.57604844 0.81741556 0. 0. ]
[ 0. 0. 1. 0. ]
[ 0. 0. 0. 1. ]]
[5.4649857 0.36596619] [[-0.40455358 -0.9145143 ]
[-0.9145143 0.40455358]]

Now as per our requirement, we need U transpose, V, and Sigma inverse, which we can compute from existing matrices. We only need to compute Sigma inverse because U and V can be directly transformed using the.T attribute. Firstly, let's convert the Sigma vector to a 2d matrix.

Input:# Transforming vector to 2d matrix
Sigma = np.diag(S)
Sigma
Output:array([[5.4649857 , 0. ],
[0. , 0.36596619]])

Next, we will estimate the inverse of the Sigma matrix.

Input:# Estimate sigma inverse from Sigma matrix 
Sigma_inverse = np.linalg.inv(Sigma)
Sigma_inverse
Output:array([[0.1829831 , 0. ],
[0. , 2.73249285]])

We will append two zero columns to match the row dimension of the input matrix A.

Input:# Adding zeros to match the row dimension of U hat transpose
# in the moore-penrose equation to get the pseudo-inverse
Sigma_plus = np.concatenate((Sigma_inverse, np.array([[0, 0], [0, 0]]).T), axis=1)
Sigma_plus
Output:array([[0.1829831 , 0. , 0. , 0. ],
[0. , 2.73249285, 0. , 0. ]])

Finally, we will place all the pieces together to get our pseudo-inverse as shown below:

Input:# Estimating pseudo-inverse using all the pieces
A_plus = np.dot(a=np.dot(a=VT.T, b=Sigma_plus), b=U.T)
A_plus
Output:array([[ 1.5, -2. , 0. , 0. ],
[-0.5, 1. , 0. , 0. ]])

Easy, right? But how can we use pseudo-inverse to estimate the least squares? Interesting, right?

In the next story, I will walk you through the unraveling of Linear Regression (aka Least Squares) implemented in the Sklearn library.

Final Thoughts and Closing Comments

There are some vital points many people fail to understand while they pursue their Data Science or AI journey. If you are one of them and looking for a way to counterbalance these cons, check out certification programs provided by INSAID on their website. If you liked this article, I would recommend going with the Global Certificate in Data Science because this one will cover your foundations plus machine learning algorithms (basic to advance).

& That’s it. I hope you liked the explanation of estimating Moore-Penrose Inverse and learned something valuable.

Follow me for more forthcoming articles based on Python, R, Data Science, Machine Learning, and Artificial Intelligence.

If you find this read helpful, then hit the Clap👏. Your encouragement will catalyze inspiration to keep me going and develop more valuable content.

--

--

Mukesh Kumar
Accredian

Data Scientist, having a robust math background, skilled in predictive modeling, data processing, and mining strategies to solve challenging business problems.