The Moore-Penrose Pseudoinverse¶

We saw that not all matrices have an inverse. It is unfortunate because the inverse is used to solve system of equations. In some cases, a system of equation has no solution, and thus the inverse doesn’t exist. However it can be useful to find a value that is almost a solution (in term of minimizing the error). We will see for instance how we can find the best-fit line of a set of data points with the pseudoinverse.

The Moore-Penrose Pseudoinverse¶

The Moore-Penrose pseudoinverse is a direct application of the SVD (see 2.8). But before all, we have to remind that systems of equations can be expressed under the matrix form.

As we have seen in 2.3, the inverse of a matrix $\boldsymbol{A}$ can be used to solve the equation $\boldsymbol{Ax}=\boldsymbol{b}$:

But in the case where the set of equations have 0 or many solutions the inverse cannot be found and the equation cannot be solved. The pseudoinverse is $\boldsymbol{A}^+$ such as:

minimizing

The following formula can be used to find the pseudoinverse:

with $\boldsymbol{U}$, $\boldsymbol{D}$ and $\boldsymbol{V}$ respectively the left singular vectors, the singular values and the right singular vectors of $\boldsymbol{A}$ (see the SVD in 2.8). $\boldsymbol{A}^+$ is the pseudoinverse of $\boldsymbol{A}$ and $\boldsymbol{D}^+$ the pseudoinverse of $\boldsymbol{D}$. We saw that $\boldsymbol{D}$ is a diagonal matrix and thus $\boldsymbol{D}^+$ can be calculated by taking the reciprocal of the non zero values of $\boldsymbol{D}$.

This is a bit crude but we will see some examples to clarify all of this.

Example 1.¶

Let’s see how to implement that. We will create a non square matrix $\boldsymbol{A}$, calculate its singular value decomposition and its pseudoinverse.

We can now check with the pinv() function from Numpy that the pseudoinverse is correct:

It looks good! We can now check that it is really the near inverse of $\boldsymbol{A}$. Since we know that

with

This is not bad! This is almost the identity matrix!

A difference with the real inverse is that $\boldsymbol{A}^+\boldsymbol{A}\approx\boldsymbol{I}$ but $\boldsymbol{A}\boldsymbol{A}^+\neq\boldsymbol{I}$.

Another way of computing the pseudoinverse is to use this formula:

The result is less acurate than the SVD method and Numpy pinv() uses the SVD (cf Numpy doc). Here is an example from the same matrix $\boldsymbol{A}$:

In this case the result is the same as with the SVD way.

Using the pseudoinverse to solve a overdetermined system of linear equations¶

In general there is no solution to overdetermined systems (see 2.4 ; Overdetermined systems). In the following picture, there is no point at the intersection of the three lines corresponding to three equations:

The pseudoinverse solve the system in the least square error perspective: it finds the solution that minimize the error. We will see this more explicitly with an example.

Example 2.¶

For this example we will consider this set of three equations with two unknowns:

Let’s see their graphical representation:

We actually see that there is no solution.

Putting this into the matrix form we have:

and

So we have:

We will now calculate the pseudoinverse of $\boldsymbol{A}$:

Now that we have calculated the pseudoinverse of $\boldsymbol{A}$:

we can use it to find $\boldsymbol{x}$ knowing that:

with:

So we have

In our two dimensions, the coordinates of $\boldsymbol{x}$ are

Let’s plot this point along with the equations lines:

Maybe you would have expected the point being at the barycenter of the triangle (cf. Least square solution in the triangle center). This is not the case becase the equations are not scaled the same way. Actually the point is at the intersection of the three symmedians of the triangle.

Example 3.¶

This method can also be used to fit a line to a set of points. Let’s take the following data points:

We have this set of $\boldsymbol{x}$ and $\boldsymbol{y}$ and we are looking for the line $y=mx+b$ that minimizes the error. The error can be evaluated as the sum of the differences between the fit and the actual data points. We can represent the data points with a matrix equations:

Note that here the matrix $\boldsymbol{A}$ represents the values of the coefficients. The column of 1 correspond to the intercepts (without it the fit would have the constraint to cross the origin). It gives the following set of equations:

We have the set of equations $mx+b=y$. The ones are used to give back the intercept parameter. For instance, in the first equation corresponding to the first point we have well $x=0$ and $y=2$. This can be confusing because here the vector $\boldsymbol{x}$ corresponds to the coefficients. This is because the problem is different from the other examples: we are looking for the coefficients of a line and not for $x$ and $y$ unknowns. We kept this notation to indicate the similarity with the last examples.

So we will construct these matrices and try to use the pseudoinverse to find the equation of the line minimizing the error (difference between the line and the actual data points).

Let’s start with the creation of the matrix of $\boldsymbol{A}$ and $\boldsymbol{b}$:

We can now calculate the pseudoinverse of $\boldsymbol{A}$:

and apply it to the result to find the coefficients with the formula:

These are the parameters of the fit. The slope is $m=0.21538462$ and the intercept is $b=2.2$. We will plot the data points and the regression line:

If you are not sure about the result. Just check it with another method. For instance, I double-checked with R:

a <- data.frame(x=c(0, 1, 2, 3, 3, 4),
                y=c(2, 4, 0, 2, 5, 3))

ggplot(data=a, aes(x=x, y=y)) +
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  xlim(-1, 5) +
  ylim(-1, 6)

outputs:

You can also do the fit with the Numpy polyfit() to check the parameters:

That’s good! We have seen how to use the pseudoinverse in order to solve a simple regression problem. Let’s see with a more realistic case.

Example 4.¶

To see the process with more data points we can generate data (see this nice blog post for other methods of fitting).

We will generate a column vector (see reshape() bellow) containing 100 points with random $x$ values and pseudo-random $y$ values. The function seed() from the Numpy.random package is used to freeze the randomisation and be able to reproduce the results:

We will create the matrix $\boldsymbol{A}$ from $\boldsymbol{x}$ by adding a column of ones exactly like we did in the example 3.

We can now find the pseudoinverse of $\boldsymbol{A}$ and calculate the coefficients of the regression line:

We can finally draw the point and the regression line:

Looks good!

Conclusion¶

You can see that the pseudoinverse can be very useful for this kind of problems! The series is not completely finished since we still have 3 chapters to cover. However, we have done the hardest part! We will now see two very light chapters before going to a nice example using all the linear algebra we have learn: the PCA.

References¶

Intuition¶

• Sean Owen - Pseudoinverse intuition

Numpy¶

• Numpy - linalg.pinv

• Numpy random seed

Systems of equations¶

• Using the Moore-Penrose Pseudoinverse to Solve Linear Equations

• Overdetermined systems

• Least square solution in the triangle center

• Symmedian

Least square fit¶

• Least square fitting