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
- positive definite (i.e. symmetric and only has positive eigenvalues, or symmetric and only has positive pivots)
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
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.