Pseudoinverse for Least Squares

1. What is a Pseudoinverse?

We know that to solve a simple equation like $ax=b$, we multiply by the inverse: $x=a^{-1}b$. Similarly, for a square and invertible matrix $A$, the solution to $A\mathbf{x}=\mathbf{b}$ is $\mathbf{x}=A^{-1}\mathbf{b}$.

But what happens when $A$ is not square or not invertible, as in our typical least squares problems? The standard inverse $A^{-1}$ doesn't exist. This is where the Moore-Penrose pseudoinverse (often just called the pseudoinverse), denoted $A^{+}$ or $A^{†}$, comes in. It's a generalization of the inverse that exists for any matrix.

The pseudoinverse always finds the "best possible" solution. What "best" means depends on the system:

• For an overdetermined system (more equations than unknowns, like in data fitting), it gives the least squares solution that minimizes $\Vert A\mathbf{x}-\mathbf{b}{\Vert }^2$.
• For an underdetermined system (fewer equations than unknowns), there are infinite solutions. The pseudoinverse gives the unique solution $\mathbf{x}$ that has the minimum Euclidean norm $\Vert \mathbf{x}\Vert$.

The following video provides an excellent introduction to why we need to generalize the concept of an inverse to handle the non-square matrices that arise in data analysis.

The key idea is that the pseudoinverse provides a unified way to handle all these cases. If a matrix is invertible, its pseudoinverse is simply its regular inverse: $A^{+}=A^{-1}$.

2. Constructing the Pseudoinverse with SVD

The Singular Value Decomposition provides a direct and numerically stable way to compute the pseudoinverse of any matrix.

Let's start with the SVD of an $m\times n$ matrix $A$:
$A=U\Sigma V^T$
If $A$ were square and invertible, we know its inverse would be $A^{-1}=V{\Sigma }^{-1}{U}^T$. The formula for the pseudoinverse looks remarkably similar:
$A^{+}=V{\Sigma }^{+}{U}^T$
The only new element is ${\Sigma }^{+}$, the pseudoinverse of the singular value matrix $\Sigma$.

So, how do we find ${\Sigma }^{+}$? The process is simple:

1. Take the transpose of $\Sigma$. If $\Sigma$ is $m\times n$, then ${\Sigma }^T$ is $n\times m$. This ensures the dimensions in the final multiplication for $A^{+}$ are correct.
2. On the diagonal of this transposed matrix, take the reciprocal of all the non-zero singular values. Leave any zeros as they are.

For example, if for a $3\times 2$ matrix $A$, the matrix $\Sigma$ is:
$\Sigma =\left[\begin{array}{cc}{\sigma }_1& 0\\ 0& {\sigma }_2\\ 0& 0\end{array}\right]$
Then ${\Sigma }^{+}$ is found by transposing and taking reciprocals of non-zero elements:
${\Sigma }^{+}=\left[\begin{array}{ccc}1/{\sigma }_1& 0& 0\\ 0& 1/{\sigma }_2& 0\end{array}\right]$
This procedure is the heart of computing the pseudoinverse. The fact that we only invert non-zero singular values is what makes the method so robust. In practice, this means we can set a small threshold and treat any singular values below it as zero, effectively filtering out numerical noise—a crucial feature for any real-world computation.

The following article gives a very clear, step-by-step guide to this construction.

3. Finding the Least Squares Solution

Once we have the pseudoinverse $A^{+}$, finding the least squares solution vector $\hat{\mathbf{x}}$ is astonishingly simple. We just multiply it by $\mathbf{b}$:
$\hat{\mathbf{x}}=A^{+}\mathbf{b}=V{\Sigma }^{+}{U}^T\mathbf{b}$
This single line replaces the two-step processes of the normal equations and QR decomposition.

The Geometric Intuition: Why It Works

This is elegant, but why does it give the least squares solution? The answer lies in the geometry of projection, which we've revisited several times. Remember, the least squares solution $\hat{\mathbf{x}}$ is the vector such that $A\hat{\mathbf{x}}$ is the orthogonal projection of $\mathbf{b}$ onto the column space of $A$.

Let's see what happens when we compute $A\hat{\mathbf{x}}$ using our new formula for $\hat{\mathbf{x}}$:
$A\hat{\mathbf{x}}=A(A^{+}\mathbf{b})=(U\Sigma V^T)(V{\Sigma }^{+}{U}^T\mathbf{b})$
Since $V$ is an orthogonal matrix, $V^{T}V=I$. The expression simplifies to:
$A\hat{\mathbf{x}}=U(\Sigma {\Sigma }^{+})U^T\mathbf{b}$
Let's look at the product $\Sigma {\Sigma }^{+}$. Using our example from before, where $\Sigma$ was $3\times 2$:
$\Sigma {\Sigma }^{+}=\left[\begin{array}{cc}{\sigma }_1& 0\\ 0& {\sigma }_2\\ 0& 0\end{array}\right]\left[\begin{array}{ccc}1/{\sigma }_1& 0& 0\\ 0& 1/{\sigma }_2& 0\end{array}\right]=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right]$
This resulting matrix is a special kind of identity matrix. It has 1s corresponding to the non-zero singular values and 0s elsewhere. If we let $r$ be the rank of $A$ (the number of non-zero singular values), then multiplying by $U(\Sigma {\Sigma }^{+})U^T$ is equivalent to multiplying by $U_r{U}_r^T$, where $U_r$ is the matrix containing the first $r$ columns of $U$.

And what is $U_r{U}_r^T$? It's the projection matrix onto the column space of A! (Since the first $r$ columns of $U$ form an orthonormal basis for the column space of A).

So, the equation $A\hat{\mathbf{x}}=(U{U}^T)\mathbf{b}$ tells us that our solution $\hat{\mathbf{x}}$ is precisely the vector that makes $A\hat{\mathbf{x}}$ equal to the projection of $\mathbf{b}$ onto the column space of $A$. This is the definition of the least squares solution.

The following video walks through this exact derivation, providing a clear explanation of how the $U{U}^T$ term acts as a projection.

This SVD-based method is the preferred approach in many modern numerical libraries because it's not only general but also the most numerically stable and robust of the methods we've seen.

Conclusion

In this lesson, we've brought together several key ideas from the course—SVD, orthogonality, and projections—to develop a powerful and elegant method for solving least squares problems.

Key Takeaways:

• The Moore-Penrose pseudoinverse $A^{+}$ generalizes the concept of a matrix inverse to any $m\times n$ matrix.
• It provides the least squares solution for overdetermined systems and the minimum norm solution for underdetermined systems.
• The pseudoinverse is computed robustly from the SVD via the formula $A^{+}=V{\Sigma }^{+}{U}^T$, where ${\Sigma }^{+}$ is formed by transposing $\Sigma$ and inverting its non-zero diagonal entries.
• The solution $\hat{\mathbf{x}}=A^{+}\mathbf{b}$ works because it finds the vector that causes $A\hat{\mathbf{x}}$ to be the orthogonal projection of $\mathbf{b}$ onto the column space of $A$.
• Compared to the normal equations, this method is far more numerically stable and is preferred in practice.

Preview of the Next Lesson:
We have now completed our module on applications in data analysis and seen how central linear algebra is to machine learning. Our next module, "Tensor Products for Composite Systems," will pivot towards your goal of understanding quantum mechanics. We will learn about a new way to combine vector spaces called the tensor product, which is fundamental for describing systems composed of multiple parts, like the multi-qubit systems in a quantum computer.