Three simple matrix decompositions in R

In the confusing world of matrix decompositions, three of the most useful are $A = LU$, $A =  U^T U$ (the Cholesky decomposition) and $A = QR$.

Here’s some R code to explore these useful beasts.

1. A = LU Link to heading

Present in every textbook in linear algebra, performing LU decompositions in R is quick and painless.

The easiest way to do this is to use the lu function in the Matrix package. This will give you a slight variant on the $A = LU$ decomposition. You’ll get $A = PLU$ instead, where $P$ is a permutation matrix.

Note LU decomposition can only be performed on invertible matrices.

library(Matrix)
mat = matrix(c(1,2,3,4,5,6,5,4,3,2,1,2,-3,2,5,20), ncol = 4)
det(mat) # -144 => matrix is not singular 
mat_lu = lu(mat)
expand(mat_lu) #contains P, L, U

mat2 = matrix(1,nrow = 3,ncol = 3 ) #example of a non-invertible matrix
lu(mat2) #errors out

2.  Cholesky decomposition Link to heading

Cholesky decompositions ($A =  U^T U$) only work with matrices that are

If you’re getting a strange error when calculating the Cholesky decomposition of a matrix, it’s probably because it has negative eigenvalues.

library(Matrix)
x = matrix(c(8,5,5,4), nrow = 2)
y = matrix(c(8,6,6,4), nrow = 2)
eigen(x) # positive eigenvalues
eigen(y) # negative eigenvalues
chol(x) # returns stuff 
chol(y) # errors out

You can also find the inverse of a matrix using the Cholesky decomposition:

@Two ways to find the inverse 
x = matrix(c(8,5,5,4), nrow = 2)
x_chol = chol(x)

#usual way of finding an inverse 
solve(x)
#using the Cholesky decomposition 
chol2inv(x_chol) 

#check the two are the same 
solve(x) - chol2inv(x_chol) #pretty much 0

One use of this is to generate correlated data:

#correlated random variables 
R = matrix(cbind(1,.40,.60,  .40,1,.90,  .25,.90,1),nrow=3) #symmetric matrix of correlations
U = t(chol(R))
N = matrix(rnorm(nrow(U)*100), nrow = 3, ncol = 100) #100 data points
X = U %*% N
X = as.data.frame(t(X))
library(GGally)
ggpairs(X) #check out that correlated data

3. A = QR decomposition Link to heading

This decomposition factors $A$ into matrices $Q$ and $R$, where $Q$ is an orthogonal matrix and $R$ is an upper triangular matrix.

mat = toeplitz(c(2,5,1,2,4)) #construct a matrix 
mat_qr = qr(mat)

#find Q 
qr.Q(mat_qr)
#proving that Q is orthogonal by showing QQ^T = I 
qr.Q(mat_qr)%*%t(qr.Q(mat_qr))  # is indeed the identity matrix 

#find R 
qr.R(mat_qr) #upper triangular

QR decompositions can be used to find eigenvalues with the QR algorithm.