Matrix Decompositions

Many complex matrix operations cannot be solved efficiently or with stability using the limited precision of computers. 

Matrix decompositions are methods that reduce a matrix into constituent parts that make it easier to calculate more complex matrix operations. 

Matrix decomposition methods, also called matrix factorization methods, are a foundation of linear algebra in computers, even for basic operations such as solving systems of linear equations, calculating the inverse, and calculating the determinant of a matrix.

This tutorial is divided into 4 parts; they are:

• What is a Matrix Decomposition
• LU Decomposition
• QR Decomposition
• Cholesky Decomposition

A. What is a Matrix Decomposition

A matrix decomposition is a way of reducing a matrix into its constituent parts.

It is an approach that can simplify more complex matrix operations that can be performed on the decomposed matrix rather than on the original matrix itself. 

A common analogy for matrix decomposition is the factoring of numbers, such as the factoring of 10 into 2 × 5. For this reason, matrix decomposition is also called matrix factorization.

Like factoring real values, there are many ways to decompose a matrix, hence there are a range of different matrix decomposition techniques.

Two simple and widely used matrix decomposition methods are the LU matrix decomposition and the QR matrix decomposition. 

B. LU Decomposition

The LU decomposition is for square matrices and decomposes a matrix into L and U components. A = L.U.

Where A is the square matrix that we wish to decompose, L is the lower triangle matrix and U is the upper triangle matrix.

The LU decomposition is found using an iterative numerical process and can fail for those matrices that cannot be decomposed or decomposed easily. 

A variation of this decomposition that is numerically more stable to solve in practice is called the LUP decomposition, or the LU decomposition with partial pivoting.  A = L.U.P.

The rows of the parent matrix are re-ordered to simplify the decomposition process and the additional P matrix specifies a way to permute the result or return the result to the original order.

There are also other variations of the LU. The LU decomposition is often used to simplify the solving of systems of linear equations, such as finding the coefficients in a linear regression, as well as in calculating the determinant and inverse of a matrix.

# LU decomposition
from numpy import array
from scipy.linalg import lu
# define a square matrix
A = array([
[1, 2, 3],
[4, 5, 6],
[7, 8, 9]])
print(A)
# factorize

P, L, U = lu(A)
print(P)
print(L)
print(U)
# reconstruct
B = P.dot(L).dot(U)
print(B)

-----Result-----

[[1 2 3]

[4 5 6]

[7 8 9]]

[[ 0. 1. 0.]

[ 0. 0. 1.]

[ 1. 0. 0.]]

[[ 1.        0.   0. ]

[ 0.14285714 1.   0. ]

[ 0.57142857 0.5  1. ]]

[[ 7.00000000e+00 8.00000000e+00 9.00000000e+00]

[ 0.00000000e+00 8.57142857e-01 1.71428571e+00]

[ 0.00000000e+00 0.00000000e+00 -1.58603289e-16]]

[[ 1. 2. 3.]

[ 4. 5. 6.]

[ 7. 8. 9.]]

C. QR Decomposition

The QR decomposition is for n × m matrices (not limited to square matrices) and decomposes a matrix into Q and R components. A = Q · R

Where A is the matrix that we wish to decompose, Q a matrix with the size m×m, and R is an upper triangle matrix with the size m × n. 

The QR decomposition is found using an iterative numerical method that can fail for those matrices that cannot be decomposed, or decomposed easily. 

Like the LU decomposition, the QR decomposition is often used to solve systems of linear equations, although is not limited to square matrices.

# QR decomposition
from numpy import array
from numpy.linalg import qr
# define rectangular matrix
A = array([
[1, 2],
[3, 4],
[5, 6]])
print(A)
# factorize
Q, R = qr(A, 'complete')
print(Q)
print(R)
# reconstruct
B = Q.dot(R)
print(B)

-----Result-----

[[1 2]

[3 4]

[5 6]]

[[-0.16903085  0.89708523  0.40824829]

[-0.50709255   0.27602622 -0.81649658]

[-0.84515425  -0.34503278  0.40824829]]

[[-5.91607978 -7.43735744]

[ 0.           0.82807867]

[ 0.           0. ]]

[[ 1. 2.]

[ 3. 4.]

[ 5. 6.]]

C. Cholesky Decomposition

The Cholesky decomposition is for square symmetric matrices where all values are greater than zero, so-called positive definite matrices. 

For our interests in machine learning, we will focus on the Cholesky decomposition for real-valued matrices and ignore the cases when working with complex numbers. 

The decomposition is defined as follows: A = L · LT

Where A is the matrix being decomposed, L is the lower triangular matrix and LT is the transpose of L. 

The decompose can also be written as the product of the upper triangular matrix, for example: A = UT · U

Where U is the upper triangular matrix. The Cholesky decomposition is used for solving linear least squares for linear regression, as well as simulation and optimization methods. 

When decomposing symmetric matrices, the Cholesky decomposition is nearly twice as efficient as the LU decomposition and should be preferred in these cases.

# Cholesky decomposition
from numpy import array
from numpy.linalg import cholesky
# define symmetrical matrix
A = array([
[2, 1, 1],
[1, 2, 1],
[1, 1, 2]])
print(A)
# factorize
L = cholesky(A)
print(L)
# reconstruct
B = L.dot(L.T)
print(B)

-----Result-----

[[2 1 1]

[1 2 1]

[1 1 2]]

[[ 1.41421356 0. 0. ]

[ 0.70710678 1.22474487 0. ]

[ 0.70710678 0.40824829 1.15470054]]

[[ 2. 1. 1.]

[ 1. 2. 1.]

[ 1. 1. 2.]]