Matrix Factorizations

Factoring is always fun and sometimes useful. We factor numbers and polynomials:

$$\large 312 = 2^3\cdot 3 \cdot 13 \quad \quad x^5-5x^3+4x=x(x-1)(x+1)(x-2)(x+2)$$

Turns out, we can sometimes factor matrices as well. There are many matrix factorizations serving various purposes. In 2000, the Institute of Electrical and Electronics Engineers (IEEE) made a list of the top 10 algorithms of the 20th century. One of the top 10 is the general idea of Matrix Factorization because of the many applications and advances due to these factorizations.

We will introduce three matrix factorizations today. The algorithm to compute each is different, but they are all built-in functions in Python, which is where we'll turn to actually use them.

LU Factorization

Many square matrices can be factored into the product of a lower triangular matrix times an upper triangular matrix: $\large A=LU$. For example:

$$\large \begin{bmatrix} 2 & 1 & -1 & 3 \\ -2 & 0 & 0 & 0 \\ 4 & 1 & -2 & 6 \\ -6 & -1 & 2 & -3 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 & 0 \\ -1 & 1 & 0 & 0 \\ 2 & -1 & 1 & 0 \\ -3 & 2 & -1 & 1 \end{bmatrix} \begin{bmatrix} 2 & 1 & -1 & 3 \\ 0 & 1 & -1 & 3 \\ 0 & 0 & -1 & 3 \\ 0 & 0 & 0 & 3 \end{bmatrix} $$

$\large LU$ factorization is primarily used to speed up the process of using Gaussian elimination to solve $\large A\vec{x}=\vec{b}$. If we have the $\large LU$ factorization of $\large A$, then we can write $\large A\vec{x}=\vec{b}$ as:

$$\large LU\vec{x}=\vec{b}$$

We can quickly solve this system in two steps:

1. Solve $\large L\vec{y}=\vec{b}$ for $\large \vec{y}$. The lower triangular structure makes this quick and easy.
2. Solve $\large U\vec{x}=\vec{y}$ for our solution $\large \vec{x}$.

Here's an example. Let's solve

$$\large \begin{bmatrix} 3 & -1 & 2 \\ 6 & -1 & 1 \\ -9 & 5 & -11 \end{bmatrix}\begin{bmatrix} x \\ y \\ z \\ \end{bmatrix} =\begin{bmatrix} -5 \\ 1 \\ 34 \\ \end{bmatrix} $$

Using the factorization

$$\large \begin{bmatrix} 3 & -1 & 2 \\ 6 & -1 & 1 \\ -9 & 5 & -11 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ -3 & 2 & 1 \end{bmatrix} \begin{bmatrix} 3 & -1 & 2 \\ 0 & 1 & -3 \\ 0 & 0 & 1 \end{bmatrix} $$

The process is illustrated in the video:

QR Factorization

To find the eigenvalues of an $\large n\times n$ matrix, we must find the roots of an $\large n^{th}$ degree polynomial. This is not hard if $\large n$ is small, but it's a mess if $\large n$ is large! Methods for finding roots of polynomials are slow and inaccurate.

So how does Python quickly give us the eigenvalues and eigenvectors for large matrices? Every square matrix $\large A$ can be factored into the product of an orthogonal matrix $\large Q$ and an upper triangular matrix $\large R$. Here is an example of a $\large QR$ factorization:

$$\large A=QR \quad \rightarrow \quad \begin{bmatrix} 1 & 1 & 1 & 1 \\ 1 & 2 & 1 & 1 \\ 1 & 1 & 0 & 1 \\ 1 & 2 & 2 & 3 \end{bmatrix} = \begin{bmatrix} \frac{1}{2} & -\frac{1}{2} & \frac{1}{2} & -\frac{1}{2} \\ \frac{1}{2} & \frac{1}{2} & -\frac{1}{2} & -\frac{1}{2} \\ \frac{1}{2} & -\frac{1}{2} & -\frac{1}{2} & \frac{1}{2} \\ \frac{1}{2} & \frac{1}{2} & \frac{1}{2} & \frac{1}{2} \end{bmatrix} \begin{bmatrix} 2 & 3 & 2 & 3 \\ 0 & 1 & 1 & 1 \\ 0 & 0 & 1 & 1 \\ 0 & 0 & 0 & 1 \end{bmatrix} $$

Note

The $\large QR$ factorization of a matrix is found using the Gram-Schmidt process (to generate the orthogonal columns of $\large Q$). We won't do this by hand!

What does this have to do with eigenvalues? With the $\large A=QR$ factorization, we can find the eigenvalues of $\large A$ in the following way:

1. Start with $\large A_0=A$.
2. In a loop, factor $\large A_k=Q_kR_k$.
3. Define our next $\large A$ by reversing the order of multiplication: $\large A_{k+1}=R_kQ_k$.

If we repeat this enough, the $\large A_k$ matrices will approach an upper triangular matrix with the eigenvalues of $\large A$ along the diagonal.

Note

There is a mysterious "reverse the order to RQ" step in the middle of this process. This is actually something called a "similarity transformation" because $\large A_{k+1}$ really equals $\large Q_k^{-1}A_kQ_k$, and it results in a matrix $\large A_{k+1}$ that has the same eigenvalues as $\large A_k$. So, as we generate new $\large A_k$'s, we keep the same eigenvalues. When we finally get to an upper triangular matrix, we can just read the eigenvalues off the diagonal. It's clever!

Below, we have a matrix $\large A$ as well as $\large A_5$, and $\large A_{10}$. We can see the upper triangular matrix forming, and we have our approximation of the eigenvalues (-2,-2, and 1):

$$\large A=\begin{bmatrix} 1 & 3 & 3 \\ -3 & -5 & -3\\ 3 & 3 & 1\\ \end{bmatrix}\quad \rightarrow \quad A_5=\begin{bmatrix} -1.91 & -0.06 & -7.46 \\ 0.05 & -2.03 & -4.09\\ -0.04 & 0.02 & 0.94\\ \end{bmatrix}\quad \rightarrow \quad A_{10}=\begin{bmatrix} -2.00 & 0.00 & 7.34 \\ 0.00 & -2.00 & 4.25\\ 0.00 & 0.00 & 1.00\\ \end{bmatrix} $$

We can get as close to the actual eigenvalues as we like by iterating sufficiently far in this process!

Diagonalization

On occasion, we want to compute $\large A^k$ for large $\large k$. This can be a computationally intensive process in general, but it's quite easy if $\large A$ is diagonal:

$$\large \begin{bmatrix} 2 & 0 & 0\\ 0 & -3 & 0 \\ 0 & 0 & -1\\ \end{bmatrix}^{10}=\begin{bmatrix} 2^{10}& 0 & 0\\ 0 & (-3)^{10} & 0 \\ 0 & 0 & (-1)^{10}\\ \end{bmatrix} $$

Our final factorization today allows us to write a matrix in terms of a diagonal matrix.

Definition: A square matrix $\large A$ is diagonalizable if $\large A$ can be factored to $\large A=PDP^{-1}$ for a diagonal matrix $\large D$ and an invertible matrix $\large P$.

Theorem: An $\large n\times n$ matrix $\large A$ is diagonalizable if and only if $\large A$ has $\large n$ linearly independent eigenvectors. Then the columns of $\large P$ are the eigenvectors of $\large A$, and the diagonal entries of $\large D$ are the corresponding eigenvalues.

As an example, find the diagonalization of $\large A$:

$$\large A=\begin{bmatrix} 7 & 2\\ -4 & 1\\ \end{bmatrix} $$

The solution is in the video:

Warning

The theorem says $\large A$ is diagonalizable if and only if it has a full complement of linearly independent eigenvectors. When does this not happen? When an eigenvalue is repeated with multiplicity $\large k$, there may or may not be $\large k$ linearly independent eigenvectors for that eigenvalue. It depends on how many free variables are in the eigenvector equation $\large \left(A-\lambda I\right)\vec{v}=\vec{0}$.

Powers of $\large A$

Now, notice that if $\large A=PDP^{-1}$, then

$$\large \begin{split} A^2 & = \left(PDP^{-1}\right)\left(PDP^{-1}\right)=PD\left(P^{-1}P\right)DP^{-1} = PDDP^{-1}=PD^2P^{-1}\\ A^3 & =\left(PDP^{-1}\right)A^2=\left(PDP^{-1}\right)\left(PD^2P^{-1}\right)=PD^3P^{-1}\\ \end{split} $$

In general, we have

$$\large A^k=PD^kP^{-1}$$

We will see an example in the next lesson where raising $\large A$ to a high power is necessary!

Other Factorizations?

There are many other factorizations that are useful in various situations. We've highlighted three of them here as examples of the entire idea. As you can see, they are used for fairly different purposes.

Section Outcomes

I can:

• Solve a small system given its $\large LU$ factorization.
• Find the approximate eigenvalues of $\large A$ from the output of a $\large QR$ algorithm.
• Find the diagonalization of a square matrix $\large A$.
• Raise $\large A$ to a high power given the diagonalization of $\large A$.