The pseudoinverse of a matrix is used when the standard inverse is not defined. This happens when a matrix is rectangular, or when it is square but singular. So what do you do in that situation if you still need to solve a system, find an approximate solution, or get a stable result in numerical computations? That’s exactly where the pseudoinverse becomes a practical tool. In this article, we will go through the definition, the key properties, and how to compute it using singular value decomposition (SVD). We will also explain where the entries of the matrices \( U \), \( \Sigma \), and \( V \) come from.
Pseudoinverse of a Matrix: When You Need It and What It Means
Let’s start with the basics. For a non-singular square matrix \( A \), the inverse \( A^{-1} \) exists, and we have:
\[
A^{-1} \cdot A = I,\qquad A \cdot A^{-1}=I.
\]
But if \( A \) is rectangular (\( m\times n \) with \( m\neq n \)), or if it is square but \( \det(A)=0 \), then the classical inverse does not exist. And that’s natural: in such cases, it is impossible to satisfy both requirements \( A \cdot A^{-1}=I \) and \( A^{-1} \cdot A=I \) at the same time.
So we introduce a generalization—a matrix that preserves the most important “inversion-like” properties in a correct form. This is how the Moore–Penrose pseudoinverse appears, and it is denoted by \( A^{+} \).
Formally, the pseudoinverse \( A^{+} \) of a matrix \( A\in\mathbb{R}^{m\times n} \) is the unique matrix \( A^{+}\in\mathbb{R}^{n\times m} \) that satisfies four conditions:
\[
A \cdot A^{+} \cdot A = A,\qquad
A^{+} \cdot A \cdot A^{+} = A^{+},\qquad
(A \cdot A^{+})^{*}=A \cdot A^{+},\qquad
(A^{+} \cdot A)^{*}=A^{+} \cdot A.
\]
Here \( ^{*} \) means conjugate transpose: in the complex case \( A^{*}=\overline{A}^{T} \), and for real matrices \( A^{*}=A^{T} \).
One more important point: if \( A \) is square and non-singular, then the pseudoinverse matches the usual inverse:
\[
A^{+}=A^{-1}.
\]
So this concept truly generalizes the familiar inverse and does not contradict it.
Pseudoinverse of a Matrix: Useful Properties for Solving Systems
Now a natural question comes up: what exactly do these four conditions give us? They “lock in” consistent behavior of \( A^{+} \) in problems where the classical inverse is not available.
First, we have the relations:
\[
A \cdot A^{+} \cdot A=A,\qquad A^{+} \cdot A \cdot A^{+}=A^{+}.
\]
They mean that after applying \( A^{+} \) and then multiplying by \( A \), we recover the original data in the correct form.
Second, the symmetry (or Hermitian property) of the products
\[
(A \cdot A^{+})^{*}=A \cdot A^{+},\qquad (A^{+} \cdot A)^{*}=A^{+} \cdot A
\]
guarantees that many geometric and numerical properties behave correctly. For real matrices, this simply means these products are symmetric.
Next comes the most practical link to least squares problems. If the system \( A \cdot x=b \) has no exact solution (which often happens for rectangular matrices), we look for a vector \( x \) that minimizes the error (usually in the Euclidean norm):
\[
\| A \cdot x-b \|\to \min.
\]
In this case, the expression \( x^{\star}=A^{+} \cdot b \) gives a solution to the least squares problem. And if there are multiple such solutions, then \( A^{+} \cdot b \) gives the one with the smallest norm \( \| x \| \). That’s very handy for checks and for real computations, isn’t it?
Pseudoinverse of a Matrix: Computing It via Singular Value Decomposition
The most common and most numerically stable way to find \( A^{+} \) is to use singular value decomposition (SVD). It exists for any matrix \( A\in\mathbb{R}^{m\times n} \) and has the form:
\[
A = U \cdot \Sigma \cdot V^{T},
\]
with the following sizes:
\[
U\in\mathbb{R}^{m\times m},\quad
\Sigma\in\mathbb{R}^{m\times n},\quad
V\in\mathbb{R}^{n\times n}.
\]
The matrices \( U \) and \( V \) are orthogonal:
\[
U^{T} \cdot U=I,\qquad V^{T} \cdot V=I.
\]
The matrix \( \Sigma \) is a rectangular “diagonal” matrix: on its main diagonal there are nonnegative numbers \( \sigma_i \), called singular values. There can be at most
\[
t=\min(m,n),\qquad \sigma_1\ge \sigma_2\ge \cdots \ge \sigma_t\ge 0,
\]
and all other entries of \( \Sigma \) are zero.
The standard way to write \( \Sigma \) depends on whether \( m \) or \( n \) is larger:
\[
\Sigma=
\begin{cases}
\begin{pmatrix}
D\\
0
\end{pmatrix}, & m\ge n,\\[12pt]
\begin{pmatrix}
D & 0
\end{pmatrix}, & m<n,
\end{cases}
\qquad
D=\operatorname{diag}(\sigma_1,\sigma_2,\ldots,\sigma_t).
\]
If \( r=\operatorname{rank}(A) \), then
\[
\sigma_1\ge \cdots \ge \sigma_r>0,\qquad \sigma_{r+1}=\cdots=\sigma_t=0.
\]
Once the SVD is computed, the pseudoinverse is found using the formula:
\[
A^{+}=V \cdot \Sigma^{+} \cdot U^{T}.
\]
Here \( \Sigma^{+}\in\mathbb{R}^{n\times m} \) is the matrix where each nonzero \( \sigma_i \) is replaced by \(\frac{1}{\sigma_i} \), while zeros stay zeros:
\[
\Sigma^{+}=
\begin{cases}
\begin{pmatrix}
D^{+} & 0
\end{pmatrix}, & m\ge n,\\[12pt]
\begin{pmatrix}
D^{+}\\
0
\end{pmatrix}, & m<n,
\end{cases}
\qquad
D^{+}=\operatorname{diag}(\sigma_1^{+},\sigma_2^{+},\ldots,\sigma_t^{+}),
\]
where
\[
\sigma_i^{+}=
\begin{cases}
\dfrac{1}{\sigma_i}, & \sigma_i>0,\\[6pt]
0, & \sigma_i=0.
\end{cases}
\]
So the main task is to correctly construct \( U \), \( \Sigma \), and \( V \). After that, the formula for \( A^{+} \) is applied directly.
Elements of the Singular Value Decomposition: Where They Come From
Alright, but where do the matrices \( U \) and \( V \) actually come from? To answer that, we look at two symmetric matrices:
\[
A \cdot A^{T}\in\mathbb{R}^{m\times m},\qquad A^{T} \cdot A\in\mathbb{R}^{n\times n}.
\]
They are positive semidefinite, so their eigenvalues \( \lambda_i \) are nonnegative, and their eigenvectors can be chosen to be orthonormal. One important detail: we take eigenvectors of \( A \cdot A^{T} \) and \( A^{T} \cdot A \), not of the matrix \( A \) itself. That’s exactly why SVD is so convenient to build through the spectral decomposition of these two symmetric matrices.
The key connection is this:
- The columns of \( U \) are orthonormal eigenvectors of the matrix \( A \cdot A^{T} \).
- The columns of \( V \) are orthonormal eigenvectors of the matrix \( A^{T} \cdot A \).
Let’s write that in formulas. If
\[
A \cdot A^{T} \cdot u_i=\lambda_i \cdot u_i,\qquad i=1,2,\ldots,m,
\]
then
\[
U=\begin{pmatrix}u_1 & u_2 & \cdots & u_m\end{pmatrix}.
\]
Similarly, if
\[
A^{T} \cdot A \cdot v_i=\lambda_i \cdot v_i,\qquad i=1,2,\ldots,n,
\]
then
\[
V=\begin{pmatrix}v_1 & v_2 & \cdots & v_n\end{pmatrix}.
\]
Now for the most important part—how singular values appear. The nonzero eigenvalues of \( A \cdot A^{T} \) and \( A^{T} \cdot A \) are the same (counting multiplicities). And the singular values are defined from them like this:
\[
\sigma_i=\sqrt{\lambda_i},\qquad \sigma_i\ge 0.
\]
So we first find \( \lambda_i \), and then take square roots.
How are the eigenvectors ui and vi related?
For each \( \sigma_i>0 \), we have:
\[
A \cdot v_i=\sigma_i \cdot u_i,\qquad A^{T} \cdot u_i=\sigma_i \cdot v_i.
\]
So if \( v_i \) and \( \sigma_i>0 \) are known, we can compute
\[
u_i=\frac{1}{\sigma_i} \cdot A \cdot v_i,
\]
and if \( u_i \) is known, then
\[
v_i=\frac{1}{\sigma_i} \cdot A^{T} \cdot u_i.
\]
This is exactly how, in practice, we build consistent matrices \( U \) and \( V \), and then form \( \Sigma \) from the singular values \( \sigma_i \).
After that, the final step looks familiar:
\[
A^{+}=V \cdot \Sigma^{+} \cdot U^{T}.
\]
So, constructing the pseudoinverse via SVD is based on the eigenvalues and eigenvectors of the two symmetric matrices \( A \cdot A^{T} \) and \( A^{T} \cdot A \), plus one simple rule for \( \Sigma^{+} \): invert the nonzero singular values and keep the zeros as zeros.
Pseudoinverse of a Matrix: Practical Part with Examples
To understand the method better, let’s move to practical examples. Step by step, we will work through a few typical problems where we apply what we’ve learned. Try to pay attention not only to the final answers, but also to the logic behind each step—this will help you feel more confident with similar exercises.
Example 1. What steps do you need to follow to find the pseudoinverse matrix?
To find the pseudoinverse matrix \( A^{+} \) using singular value decomposition, you usually follow these steps:
- For the given matrix \( A\in\mathbb{R}^{m\times n} \), compute two symmetric matrices:
\[
A \cdot A^{T}\in\mathbb{R}^{m\times m},\qquad A^{T} \cdot A\in\mathbb{R}^{n\times n}.
\]
- For the matrix \(A \cdot A^{T} \), find the eigenvalues \( \lambda_i \) and orthonormal eigenvectors \( u_i \). Build the matrix
\[
U=\begin{pmatrix}u_1 & u_2 & \cdots & u_m\end{pmatrix}.
\]
- For the matrix \(A^{T} \cdot A \), find the eigenvalues \( \lambda_i \) and orthonormal eigenvectors \( v_i \). Build the matrix
\[
V=\begin{pmatrix}v_1 & v_2 & \cdots & v_n\end{pmatrix}.
\]
- From the eigenvalues, obtain the singular values \( \sigma_i=\sqrt{\lambda_i} \) (\( \sigma_i\ge 0 \)). Then form the matrix \( \Sigma\in\mathbb{R}^{m\times n} \) so that the nonzero \( \sigma_i \) are on the main diagonal, and all other entries are zero.
- Build \( \Sigma^{+}\in\mathbb{R}^{n\times m} \): place the reciprocals of the nonzero singular values on its diagonal, and put zeros where \( \sigma_i=0 \).
- Write the pseudoinverse using the formula:
\[
A^{+}=V \cdot \Sigma^{+} \cdot U^{T}.
\]
Example 2. Find the pseudoinverse matrix for the given matrix
\[
A=
\begin{pmatrix}
2 & 1 \\
1 & 2 \\
3 & 3
\end{pmatrix}.
\]
We start by building the symmetric matrices \( A \cdot A^{T} \) and \(A^{T} \cdot A \):
\[
\begin{gathered}
A \cdot A^{T}=
\begin{pmatrix}
2 & 1\\
1 & 2\\
3 & 3
\end{pmatrix} \cdot
\begin{pmatrix}
2 & 1 & 3\\
1 & 2 & 3
\end{pmatrix}
=
\begin{pmatrix}
5 & 4 & 9\\
4 & 5 & 9\\
9 & 9 & 18
\end{pmatrix},
\\[6pt]
A^{T} \cdot A=
\begin{pmatrix}
2 & 1 & 3\\
1 & 2 & 3
\end{pmatrix} \cdot
\begin{pmatrix}
2 & 1\\
1 & 2\\
3 & 3
\end{pmatrix}
=
\begin{pmatrix}
14 & 13\\
13 & 14
\end{pmatrix}.
\end{gathered}
\]
Next, we find the spectral data. For \( A \cdot A^{T} \), we can take the eigenvalues and corresponding eigenvectors as:
\[
\lambda_1=27,\quad x_1=
\begin{pmatrix}
1\\
1\\
2
\end{pmatrix},
\qquad
\lambda_2=1,\quad x_2=
\begin{pmatrix}
-1\\
1\\
0
\end{pmatrix},
\qquad
\lambda_3=0,\quad x_3=
\begin{pmatrix}
-1\\
-1\\
1
\end{pmatrix}.
\]
We normalize them and form the orthonormal columns of \( U \). As a result, we get:
\[
U=
\begin{pmatrix}
\frac{1}{\sqrt{1^{2}+1^{2}+2^{2}}} & \frac{-1}{\sqrt{(-1)^{2}+1^{2}+0^{2}}} & \frac{-1}{\sqrt{(-1)^{2}+(-1)^{2}+1^{2}}} \\
\frac{1}{\sqrt{1^{2}+1^{2}+2^{2}}} & \frac{1}{\sqrt{(-1)^{2}+1^{2}+0^{2}}} & \frac{-1}{\sqrt{(-1)^{2}+(-1)^{2}+1^{2}}} \\
\frac{2}{\sqrt{1^{2}+1^{2}+2^{2}}} & \frac{0}{\sqrt{(-1)^{2}+1^{2}+0^{2}}} & \frac{1}{\sqrt{(-1)^{2}+(-1)^{2}+1^{2}}}
\end{pmatrix}=
\begin{pmatrix}
0.408 & -0.707 & -0.577\\
0.408 & 0.707 & -0.577\\
0.816 & 0 & 0.577
\end{pmatrix}.
\]
For \( A^{T} \cdot A \), the eigenvalues and eigenvectors can be written as:
\[
\lambda_1=27,\quad y_1=
\begin{pmatrix}
1\\
1
\end{pmatrix},
\qquad
\lambda_2=1,\quad y_2=
\begin{pmatrix}
-1\\
1
\end{pmatrix}.
\]
After normalization, we form \( V \):
\[
V=
\begin{pmatrix}
\frac{1}{\sqrt{1^{2}+1^{2}}} & \frac{-1}{\sqrt{(-1)^{2}+1^{2}}} \\
\frac{1}{\sqrt{1^{2}+1^{2}}} & \frac{1}{\sqrt{(-1)^{2}+1^{2}}}
\end{pmatrix}=
\begin{pmatrix}
0.707 & -0.707\\
0.707 & 0.707
\end{pmatrix}.
\]
Now we move to the singular values. Since \( \sigma_i=\sqrt{\lambda_i} \), we have:
\[
\sigma_1=5.196,\qquad \sigma_2=1.
\]
So the matrix \( \Sigma\in\mathbb{R}^{3\times 2} \) becomes:
\[
\Sigma=
\begin{pmatrix}
5.196 & 0\\
0 & 1\\
0 & 0
\end{pmatrix}.
\]
Next, we build \( \Sigma^{+}\in\mathbb{R}^{2\times 3} \). In place of \( 5.196 \) we put its reciprocal, and we keep \( 1 \) and the zeros unchanged:
\[
\Sigma^{+}=
\begin{pmatrix}
0.192 & 0 & 0\\
0 & 1 & 0
\end{pmatrix}.
\]
Now we apply the key formula:
\[
A^{+}=V \cdot \Sigma^{+} \cdot U^{T}.
\]
After multiplying, we obtain the pseudoinverse matrix:
\[
A^{+}=
\begin{pmatrix}
0.707 & -0.707\\
0.707 & 0.707
\end{pmatrix} \cdot
\begin{pmatrix}
0.192 & 0 & 0\\
0 & 1 & 0
\end{pmatrix} \cdot
\begin{pmatrix}
0.408 & 0.408 & 0.816\\
-0.707 & 0.707 & 0\\
-0.577 & -0.577 & 0.577
\end{pmatrix}=
\begin{pmatrix}
0.556 & -0.444 & 0.111\\
-0.444 & 0.556 & 0.111
\end{pmatrix}.
\]
Dimension check: \( A\in\mathbb{R}^{3\times 2} \), so \( A^{+}\in\mathbb{R}^{2\times 3} \), which matches the result we obtained.
Example 3. Find the pseudoinverse matrix for the given matrix
\[
A=
\begin{pmatrix}
2 & 0 & 0\\
0 & 1 & 0
\end{pmatrix}.
\]
Just like in the previous example, we start by building the symmetric matrices \( A \cdot A^{T} \) and \( A^{T} \cdot A \):
\[
\begin{gathered}
A \cdot A^{T}=
\begin{pmatrix}
2 & 0 & 0\\
0 & 1 & 0
\end{pmatrix} \cdot
\begin{pmatrix}
2 & 0\\
0 & 1\\
0 & 0
\end{pmatrix}
=
\begin{pmatrix}
4 & 0\\
0 & 1
\end{pmatrix},
\\[6pt]
A^{T} \cdot A=
\begin{pmatrix}
2 & 0\\
0 & 1\\
0 & 0
\end{pmatrix} \cdot
\begin{pmatrix}
2 & 0 & 0\\
0 & 1 & 0
\end{pmatrix}
=
\begin{pmatrix}
4 & 0 & 0\\
0 & 1 & 0\\
0 & 0 & 0
\end{pmatrix}.
\end{gathered}
\]
Now we look at the spectral data. For \( A \cdot A^{T} \), we can take:
\[
\lambda_1=4,\quad x_1=
\begin{pmatrix}
1\\
0
\end{pmatrix},
\qquad
\lambda_2=1,\quad x_2=
\begin{pmatrix}
0\\
1
\end{pmatrix}.
\]
These vectors are already orthonormal, so the matrix \( U \) is simply:
\[
U=
\begin{pmatrix}
1 & 0 \\
0 & 1
\end{pmatrix}.
\]
For \( A^{T} \cdot A \), we have:
\[
\lambda_1=4,\quad y_1=
\begin{pmatrix}
1\\
0\\
0
\end{pmatrix},
\qquad
\lambda_2=1,\quad y_2=
\begin{pmatrix}
0\\
1\\
0
\end{pmatrix},
\qquad
\lambda_3=0,\quad y_3=
\begin{pmatrix}
0\\
0\\
1
\end{pmatrix}.
\]
These vectors are also already orthonormal, so
\[
V=
\begin{pmatrix}
1 & 0 & 0\\
0 & 1 & 0\\
0 & 0 & 1
\end{pmatrix}.
\]
Since the matrix \( A \) has size \( 2\times 3 \), the matrix \( \Sigma \) can have only \( t=\min(2,3)=2 \) singular values. So here we take:
\[
\sigma_1=2,\qquad \sigma_2=1.
\]
Then \( \Sigma\in\mathbb{R}^{2\times 3} \) is:
\[
\Sigma=
\begin{pmatrix}
2 & 0 & 0 \\
0 & 1 & 0
\end{pmatrix}.
\]
Next, we build \( \Sigma^{+}\in\mathbb{R}^{3\times 2} \). We replace the nonzero singular values with their reciprocals, and keep the zero positions as zeros:
\[
\Sigma^{+}=
\begin{pmatrix}
0.5 & 0\\
0 & 1\\
0 & 0
\end{pmatrix}.
\]
Now we apply the formula:
\[
A^{+}=V \cdot \Sigma^{+} \cdot U^{T}.
\]
Since \( V \) and \( U \) in this example are identity matrices, the multiplication does not change \( \Sigma^{+} \). Therefore:
\[
A^{+}=
\begin{pmatrix}
0.5 & 0\\
0 & 1\\
0 & 0
\end{pmatrix}.
\]
Pseudoinverse of a Matrix in Code: Try Building Your Own Mini Calculator
Imagine that the flowchart below is not just a picture, but a clear plan for a small programming project. Why not use it as a hint and write a simple program in your favorite programming language that computes the pseudoinverse of a matrix using singular value decomposition? This kind of task is great for training your attention to detail, and at the same time it shows how a mathematical idea turns into a practical tool—something you can easily test on your own examples.
And honestly, it’s satisfying to see the result not only on paper, but also as code that gives you an answer on the very first run.
