Three simple matrix decompositions in R
11 Dec 2016 · 403 words

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

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

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

• real
• symmetric
• positive definite (i.e. symmetric and only has positive eigenvalues, or symmetric and only has positive pivots)
• square

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

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.