### create a pxp Sigma
p = 5
rho = .8 
Sigma = matrix(rep(rho,p^2),ncol=p)
diag(Sigma) = 1.0


### get the cholesky
L = t(chol(Sigma))

### check Sigma = LL'
L %*% t(L)

### get the eigen values/vectors
eg = eigen(Sigma)
names(eg)

### Sigma = PDP'
D = diag(eg$values)
P = eg$vectors

### check Sigma = PDP'
check1 = Sigma %*% P
check2 = P %*% D
check1 - check2

pdpt = P %*% D %*% t(P)
Sigma - pdpt

### solve a system of equations with a lower triangular
b = matrix(1:p,ncol=1)
x = solve(L) %*% b

xx = forwardsolve(L,b)

### compare the time to invert with the time to choleski
p = 500
rho = .8
Sigma = matrix(rep(rho,p^2),ncol=p)
diag(Sigma) = 1.0

t1 = system.time({Si = solve(Sigma)})

t2 = system.time({R = chol(Sigma)})