Matrix Decomposition Demystified

Eigen Decomposition, SVD, and Pseudo-inverse Matrix Made Easy

1 Matrix Decomposition

In this article, we discuss the following Matrix Decomposition topics:

• Square Matrix
• Eigenvalue and Eigenvector
• Symmetric Matrix
• Eigen Decomposition
• Orthogonal Matrix
• Singular Value Decomposition
• PSeudo-inverse Matrix

2 Square Matrix

Eigen Decomposition works only with square matrices.

As a quick reminder, let’s look at a square matrix. In square matrices, the number of rows and columns are the same.

For example,

$$A=\left[\begin{array}{cc}4& 1\\ 3& 2\end{array}\right]$$

That was easy. Let’s continue with Eigenvalue and Eigenvector concepts.

3 Eigenvalue and Eigenvector

When a square matrix A and a vector $x$ have the following relationship:

$$A\mathbf{x}=\lambda \mathbf{x}$$

, we call λ eigenvalue and the vector x eigenvector. λ is a scalar.

The above means that multiplying vector $x$ with matrix $A$ produces the same result as multiplying vector $x$ with scalar value $\lambda$.

Let’s move everything to the left side of the equation to find the condition for having eigenvalues and eigenvectors.

$$(A-\lambda I)\mathbf{x}=\mathbf{0}$$

For vector $x$ to be non-zero, no inverse matrix of $(A—\lambda I)$ should exist. In other words, $(A — I) $’s determinant must be zero.

$$det(A-\lambda I)=0$$

Let’s use the 2x2 square matrix as an example to calculate the determinant:

$$A=\left[\begin{array}{cc}4& 1\\ 3& 2\end{array}\right]$$

So, $(A — \lambda I)$ is:

$$A-\lambda I=\left[\begin{array}{cc}4-\lambda & 1\\ 3& 2-\lambda \end{array}\right]$$

And we calculate the determinant (and we want it to be zero):

$$\begin{array}{rl}det(A-\lambda I)& =\left|\begin{array}{cc}4-\lambda & 1\\ 3& 2-\lambda \end{array}\right|\\ \\ & =(4-\lambda )(2-\lambda )-3\\ \\ & =8-6\lambda +{\lambda }^2-3\\ \\ & ={\lambda }^2-6\lambda +5\\ \\ & =(\lambda -1)(\lambda -5)\\ \\ & =0\end{array}$$

Therefore, $\lambda$ is either 1 or 5.

Let’s define the eigenvector as follows:

$$\mathbf{x}=\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]$$

As a reminder, we have the following equation for eigenvalue and eigenvector:

$$(A-\lambda I)=0$$

We can calculate the eigenvector $x$ for $\lambda =1$ as follows:

$$\begin{array}{rl}A& =\left[\begin{array}{cc}4-1& 1\\ 3& 2-1\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right]\\ \\ & \therefore 3x_1+x_2=0\end{array}$$

To satisfy the above conditions, $x_1$ and $x_2$ are:

$$\begin{array}{rl}{x}_1& =-\frac{1}{3}t\\ x_2& =t\end{array}$$

Although t can be any value, we often make the L2 norm to be a size of 1 since, otherwise, we have an infinite number of possible solutions:

$$x_1^2+x_2^2=1$$

The eigenvector for the eigenvalue $\lambda =1$ is:

$$\mathbf{x}=\left[\begin{array}{c}-\frac{1}{\sqrt{10}}\\ \frac{3}{\sqrt{10}}\end{array}\right]$$

or

$$\mathbf{x}=\left[\begin{array}{c}\frac{1}{\sqrt{10}}\\ -\frac{3}{\sqrt{10}}\end{array}\right]$$

They are the same, except that one vector direction is the opposite. So, I’ll choose the first one as the eigenvector for $\lambda =1$.

Let’s make sure this works as intended:

$$\begin{array}{rl}A\mathbf{x}& =\left[\begin{array}{cc}4& 1\\ 3& 2\end{array}\right]\left[\begin{array}{c}-\frac{1}{\sqrt{10}}\\ \frac{3}{\sqrt{10}}\end{array}\right]\\ \\ & =\left[\begin{array}{c}4\cdot -\frac{1}{\sqrt{10}}+1\cdot \frac{3}{\sqrt{10}}\\ 3\cdot -\frac{1}{\sqrt{10}}+2\cdot \frac{3}{\sqrt{10}}\end{array}\right]\\ \\ & =\left[\begin{array}{c}-\frac{1}{\sqrt{10}}\\ \frac{3}{\sqrt{10}}\end{array}\right]\\ \\ & =1\mathbf{x}\end{array}$$

We can solve for $\lambda =5$ using the same process:

$$\mathbf{x}=\left[\begin{array}{c}\frac{1}{sqrt2}\\ \frac{1}{sqrt2}\end{array}\right]$$

For 3×3 or larger matrices, it becomes too tedious to calculate all that manually. We can write a Python script with NumPy to do the same:

import numpy as np

A = np.array([
    [1, 2, 3],
    [4, 5, 6],
    [7, 8, 9]])

print('A=')
print(A)
print()

values, vectors = np.linalg.eig(A)

print('eigenvalues')
print(values)
print()

print('eigenvectors')
print(vectors)

Below is the output:

A=
[[1 2 3]
 [4 5 6]
 [7 8 9]]

eigenvalues
[ 1.61168440e+01 -1.11684397e+00 -1.30367773e-15]

eigenvectors
[[-0.23197069 -0.78583024 0.40824829]
 [-0.52532209 -0.08675134 -0.81649658]
 [-0.8186735 0.61232756 0.40824829]]

We can confirm they work as expected:

λ = values[0]
x = vectors[:, 0]

print('λ * x')
print(λ * x)
print()
print('A @ x')
print(A @ x)

Below is the output:

λ * x
[ -3.73863537 -8.46653421 -13.19443305]

A @ x
[ -3.73863537 -8.46653421 -13.19443305]

This works fine but if you want to keep the vectors in column format, you can do the following:

λ = values[0]
x = vectors[:, 0:1]

print('λ * x')
print(λ * x)
print()
print('A @ x')
print(A @ x)

Below is the output:

λ * x
[[ -3.73863537]
 [ -8.46653421]
 [-13.19443305]]

A @ x
[[ -3.73863537]
 [ -8.46653421]
 [-13.19443305]]

4 Symmetric Matrix

Eigenvectors of a symmetric matrix are orthogonal to each other. This fact will become helpful later on. So, let’s prove it here.

Suppose ${\lambda }_1$ and ${\lambda }_2$ (${\lambda }_1\ne {\lambda }_2$) are eigenvalues, and ${\mathbf{x}}_1$ and ${\mathbf{x}}_2$ are the corresponding eigenvectors:

$$\begin{array}{rl}{\lambda }_1{\mathbf{x}}_1\cdot {\mathbf{x}}_2& =({\lambda }_1{\mathbf{x}}_1{)}^{\mathrm{\top }}{\mathbf{x}}_2\\ & =(A{\mathbf{x}}_1{)}^{\mathrm{\top }}{\mathbf{x}}_2\\ & ={\mathbf{x}}_1^{\mathrm{\top }}{A}^{\mathrm{\top }}{\mathbf{x}}_2\\ & ={\mathbf{x}}_1^{\mathrm{\top }}A{\mathbf{x}}_2\text{A is symmetric}\\ & ={\mathbf{x}}_1^{\mathrm{\top }}{\lambda }_2{\mathbf{x}}_2\\ & ={\lambda }_2{\mathbf{x}}_1^{\mathrm{\top }}{\mathbf{x}}_2\\ & ={\lambda }_2{\mathbf{x}}_1\cdot {\mathbf{x}}_2\end{array}$$

Therefore,

$${\lambda }_1{\mathbf{x}}_1\cdot {\mathbf{x}}_2={\lambda }_2{\mathbf{x}}_1\cdot {\mathbf{x}}_2$$

However, ${\lambda }_1\ne {\lambda }_2$. As such, the following must be true:

$${\mathbf{x}}_1\cdot {\mathbf{x}}_2=\mathbf{0}$$

Therefore, eigenvectors of a symmetric matrix are orthogonal to each other.

Now we are ready to discuss how Eigen Decomposition works.

5 Eigen Decomposition

Using eigenvalues and eigenvectors, we can decompose a square matrix A as follows:

$$A=Q\mathrm{\Lambda }{Q}^{-1}$$

Q is a matrix that has eigenvectors in the columns.

$$Q=[{\mathbf{x}}_1{\mathbf{x}}_2\dots {\mathbf{x}}_n]$$

$\mathrm{\Lambda }$ (uppercase lambda) is a diagonal matrix that has eigenvalues as the diagonal elements:

$$\mathrm{\Lambda }=\left[\begin{array}{cccc}{\lambda }_1& 0& \dots & 0\\ 0& {\lambda }_2& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & {\lambda }_n\end{array}\right]$$

We put eigenvalues in descending order to make the diagonal matrix $\mathrm{\Lambda }$ unique:

$${\lambda }_1>{\lambda }_2>\cdots >{\lambda }_n$$

To prove that such Eigen Decomposition is possible, we adjust the equation slightly:

$$\begin{array}{rl}A& =Q\mathrm{\Lambda }{Q}^{-1}\\ AQ& =Q\mathrm{\Lambda }\end{array}$$

And we will show that $AQ=Q\mathrm{\Lambda }$ is true.

$$\begin{array}{rl}AQ& =A[{\mathbf{x}}_1{\mathbf{x}}_2\dots {\mathbf{x}}_n]\\ & =[A{\mathbf{x}}_{1}A{\mathbf{x}}_2\dots A{\mathbf{x}}_n]\\ & =[{\lambda }_1{\mathbf{x}}_1{\lambda }_2{\mathbf{x}}_2\dots {\lambda }_n{\mathbf{x}}_n]\\ & =[{\mathbf{x}}_1{\mathbf{x}}_2\dots {\mathbf{x}}_n]\left[\begin{array}{cccc}{\lambda }_1& 0& \dots & 0\\ 0& {\lambda }_2& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & {\lambda }_n\end{array}\right]\\ & =Q\mathrm{\Lambda }\end{array}$$

In other words, AQ is equivalent to each eigenvector inside matrix Q multiplied by the corresponding eigenvalue.

So, we have proved that Eigen Decomposition is possible.

Let’s try Eigen Decomposition with Numpy:

import numpy as np

A = np.array([
    [1, 2, 3],
    [4, 5, 6],
    [7, 8, 9]], dtype=np.float64)

print('A=')
print(A)
print()

# Eigenvalue and Eigenvector
values, vectors = np.linalg.eig(A)

Λ = np.diag(values)
Q = vectors

QΛQinv = Q @ Λ @ np.linalg.inv(Q)

print('QΛQinv=')
print(QΛQinv)

Below is the output:

A=
[[1. 2. 3.]
 [4. 5. 6.]
 [7. 8. 9.]]

QΛQinv=
[[1. 2. 3.]
 [4. 5. 6.]
 [7. 8. 9.]]

6 Orthogonal Matrix

As mentioned before, when matrix A is symmetric, eigenvectors in matrix Q are orthogonal to each other.

When we also make the L2 norm of all eigenvectors 1 (orthonormal), we call Q an orthogonal matrix.

When Q is an orthogonal matrix,

$$Q{Q}^{\mathrm{\top }}=I$$

As such,

$$Q^{\mathrm{\top }}=Q^{-1}$$

Therefore,

$$A=Q\mathrm{\Lambda }{Q}^{-1}=Q\mathrm{\Lambda }{Q}^{\mathrm{\top }}$$

The above fact will be helpful when we discuss the singular value decomposition.

Let’s do QΛQ^T with Numpy:

import numpy as np

# Asymmetric Matrix
A = np.array([
    [1, 2, 3],
    [2, 5, 6],
    [3, 6, 9]], dtype=np.float64)

print('A=')
print(A)
print()

# Eigenvalue and Eigenvector
values, vectors = np.linalg.eig(A)

Λ = np.diag(values)
Q = vectors

QΛQT = Q @ Λ @ Q.T

print('QΛQT=')
print(QΛQT)

Below is the output:

A=
[[1. 2. 3.]
 [2. 5. 6.]
 [3. 6. 9.]]

QΛQT=
[[1. 2. 3.]
 [2. 5. 6.]
 [3. 6. 9.]]

It’s fun to see the theory works precisely as expected.

Also, let’s check Q is an orthogonal matrix:

QQT = Q @ Q.T

print('QQT=')
print(QQT)

Below is the output:

QQT=
[[ 1.00000000e+00 5.55111512e-17 -1.80411242e-16]
 [ 5.55111512e-17 1.00000000e+00 0.00000000e+00]
 [-1.80411242e-16 0.00000000e+00 1.00000000e+00]]

Here we see the limitation of numerical computing. Non-diagonal elements should be 0; instead, some are very small values.

So far, we are dealing with square matrices as Eigen Decomposition works only with square matrices. Next, we’ll see SVD that works with non-square matrices.

7 Singular Value Decomposition

SVD exists for any matrix, even for non-square matrices.

Suppose matrix A is m×n (m ≠ n), we can still decompose matrix A as follows:

$$A=U\mathrm{\Sigma }{V}^{\mathrm{\top }}$$

• $U$ is an orthogonal matrix (m×m).
• $\mathrm{\Sigma }$ is a diagonal matrix (m×n).
• $V$ is an orthogonal matrix (n×n).

It looks too abstract — a visualization may help:

$\mathrm{\Sigma }$ has the following structure. The top-left part is a diagonal matrix with singular values (more on this later) in the diagonal elements. Other elements of $\mathrm{\Sigma }$ are all zeros.

As a convention, we sort the singular values in descending order:

$${\sigma }_1\ge {\sigma }_2\ge \dots \ge {\sigma }_s\ge 0$$

The number of singular values is the same as the rank of matrix A, which is the number of independent column vectors.

So, how can we create such a decomposition?

The main idea is that we make a square matrix like the one below:

$$A{A}^{\mathrm{\top }}$$

or

$$A^{\mathrm{\top }}A$$

We can choose either, but using the one with fewer elements is easier.

For example, if matrix A is 100×10,000. $A{A}^T$ is 100×100, and $A^{T}A$ is 10,000×10,000. So, I would choose $A{A}^T$. How about you?

Suppose we perform SVD on matrix A using $A{A}^T$:

$$\begin{array}{rl}A{A}^{\mathrm{\top }}& =U\mathrm{\Sigma }{V}^{\mathrm{\top }}(U\mathrm{\Sigma }{V}^{\mathrm{\top }}{)}^{\mathrm{\top }}\\ & =U\mathrm{\Sigma }{V}^{\mathrm{\top }}V{\mathrm{\Sigma }}^{\mathrm{\top }}{U}^{\mathrm{\top }}\text{V is an orthogonal matrix}\\ & =U\mathrm{\Sigma }{\mathrm{\Sigma }}^{\mathrm{\top }}{U}^{\mathrm{\top }}\\ & =U{\mathrm{\Sigma }}_m^2{U}^{\mathrm{\top }}\end{array}$$

where

$${\mathrm{\Sigma }}_m^2=\mathrm{\Sigma }{\mathrm{\Sigma }}^{\mathrm{\top }}$$

is a diagonal matrix (m×m) with squared sigmas as the top-left diagonal elements. These values are eigenvalues of AA^T since we have eigen-decomposed $A{A}^T$ as follows:

$$A{A}^{\mathrm{\top }}=U{\mathrm{\Sigma }}_m^2{U}^{\mathrm{\top }}$$

The column vectors in matrix U are eigenvectors because:

$$A{A}^{\mathrm{\top }}U=U{\mathrm{\Sigma }}_m^2\text{U is an orthogonal matrix}$$

In other words, if we can eigen-decompose $A{A}^T$, we can calculate U and diagonal sigmas, and then we can calculate V.

We can consider the same calculation with $A^{T}A$:

$$\begin{array}{rl}{A}^{\mathrm{\top }}A& =(U\mathrm{\Sigma }{V}^{\mathrm{\top }}{)}^{\mathrm{\top }}U\mathrm{\Sigma }{V}^{\mathrm{\top }}\\ & =V{\mathrm{\Sigma }}^{\mathrm{\top }}{U}^{\mathrm{\top }}U\mathrm{\Sigma }{V}^{\mathrm{\top }}\text{U is an orthogonal matrix}\\ & =V{\mathrm{\Sigma }}^{\mathrm{\top }}\mathrm{\Sigma }{V}^{\mathrm{\top }}\\ & =V{\mathrm{\Sigma }}_n^2{V}^{\mathrm{\top }}\end{array}$$

where

$${\mathrm{\Sigma }}_n^2={\mathrm{\Sigma }}^{\mathrm{\top }}\mathrm{\Sigma }$$

is a diagonal matrix (n×n) with squared sigmas as the top-left diagonal elements. These values are eigenvalues since we eigen-decompose $A^{T}A$ as follows:

$$A^{\mathrm{\top }}A=V{\mathrm{\Sigma }}_n^2{V}^{\mathrm{\top }}$$

The column vectors in matrix V are eigenvectors because:

$$A^{\mathrm{\top }}AV=V{\mathrm{\Sigma }}_n^2$$

In other words, if we can eigen-decompose A^TA, we can calculate V and diagonal sigmas, and then we can calculate U.

So, if we eigen-decompose either $A{A}^T$ or $A^{T}A$, we can perform singular value decomposition.

Again, we should choose a smaller one of $A{A}^T$ or $A^{T}A$ to make our life easier.

Let’s try SVD with Numpy:

import numpy as np

A = np.array([
    [1, 2],
    [3, 4],
    [5, 6]], dtype=np.float64)

print('A=')
print(A)
print()

U, S, VT = np.linalg.svd(A)

print('U=')
print(U)
print()

print('S=')
print(S)
print()

print('VT=')
print(VT)
print()

Below is the output:

A=
[[1. 2.]
 [3. 4.]
 [5. 6.]]

U=
[[-0.2298477 0.88346102 0.40824829]
 [-0.52474482 0.24078249 -0.81649658]
 [-0.81964194 -0.40189603 0.40824829]]

S=
[9.52551809 0.51430058]

VT=
[[-0.61962948 -0.78489445]
 [-0.78489445 0.61962948]]

Let’s set up matrix Σ as follows:

Σ = np.zeros_like(A, dtype=np.float64)
Σ[:2] = np.diag(S)

print('Σ=')
print(Σ)
print()

Below is the output:

Σ=
[[9.52551809 0. ]
 [0. 0.51430058]
 [0. 0. ]]

Now, we can confirm $A=U\mathrm{\Sigma }{V}^{\mathrm{\top }}$.

print('UΣVT=')
print(U @ Σ @ VT)
print()

Below is the output:

UΣVT=
[[1. 2.]
 [3. 4.]
 [5. 6.]]

Now that we know how SVD works, we are ready to work on pseudo-inverse matrices.

8 Pseudo-inverse Matrix

An inverse matrix does not always exist, even for a square matrix. However, a pseudo-inverse — also called Moore Penrose inverse — matrix exists for non-square matrices.

For example, matrix A is m×n.

$$A\mathbf{x}=\mathbf{y}$$

Using a pseudo-inverse matrix $A^{+}$, we can perform the following conversion:

$$\mathbf{x}=A^{+}\mathbf{y}$$

We define a pseudo-inverse matrix $A^{+}$ as:

$$A+=V{D}^{+}{U}^{\mathrm{\top }}$$

V and U are from SVD:

$$A=U\mathrm{\Sigma }{V}^{\mathrm{\top }}$$

We make $D^{+}$ by transposing $\mathrm{\Sigma }$ and inverse all the diagonal elements. Suppose $\mathrm{\Sigma }$ is defined as follows:

$$\mathrm{\Sigma }=\left[\begin{array}{cccccc}{\sigma }_1& & & & & \\ & {\sigma }_2& & & & \\ & & \ddots & & & & \\ & & & {\sigma }_s& & & \\ & & & & 0& & \\ & & & & & \ddots & \\ & & & & & & 0\end{array}\right]$$

Then $D^{+}$ is defined as follows:

$$D^{+}=\left[\begin{array}{cccccc}\frac{1}{{\sigma }_1}& & & & & \\ & \frac{1}{{\sigma }_2}& & & & \\ & & \ddots & & & & \\ & & & \frac{1}{{\sigma }_s}& & & \\ & & & & 0& & \\ & & & & & \ddots & \\ & & & & & & 0\end{array}\right]$$

Now, we can see how $A^{+}A$ works:

$$\begin{array}{rl}{A}^{+}A& =(V{D}^{+}{U}^{\mathrm{\top }})(U\mathrm{\Sigma }{V}^{\mathrm{\top }})\\ & =V{D}^{+}\mathrm{\Sigma }{V}^{\mathrm{\top }}\text{Note: }{D}^{+}\mathrm{\Sigma }=I\\ & =V{V}^{\mathrm{\top }}\\ & =I\end{array}$$

In the same way, $A{A}^{+}=I$.

In summary, if we can perform SVD on matrix A, we can calculate A^+ by VD^+UT, which is a pseudo-inverse matrix of A.

Let’s try Pseudo-inverse with Numpy:

import numpy as np

A = np.array([
    [1, 2],
    [3, 4],
    [5, 6]], dtype=np.float64)

AP = np.linalg.pinv(A)

print(‘AP @ A’)
print(AP @ A)

The below is the output:

AP @ A
[[ 1.0000000000000004e+00 0.0000000000000000e+00]
 [-4.4408920985006262e-16 9.9999999999999956e-01]]

Congratulations! You’ve made it this far. I hope you’ve enjoyed it!

9 References

• Eigendecomposition of a matrix
  Wikipedia
• Singular Value Decomposition
  Wikipedia
• Moore Penrose inverse
  Wikipedia
• Up-sampling with Transposed Convolution