# 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

- 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 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.