sessionInfo()

library(Matrix)
library(mvtnorm)
library(microbenchmark)

A <- t(Matrix(c(4, 12, -16, 12, 37, -43, -16, -43, 98), nrow=3))  # check out this is pd!
A

# Cholesky without pivoting
Achol <- Matrix::chol(Matrix::symmpart(A))
Achol

class(Achol)

getSlots("Cholesky")

str(Achol)

b <- c(1.0, 2.0, 3.0)
solve(A, b) # this does LU; wasteful!; 2/3 n^3 + 2n^2

y <- solve(t(Achol), b) # two triangular solves; only 2n^2 flops
solve(Achol, y)

det(A) # this actually does LU; wasteful!

det(Achol)^2 # cheap

solve(A) # this does LU!

solve(Achol, solve(t(Achol))) # 2n triangular solves

(A <- Matrix(c(1,1,.5,1,1,.5,.5,.5,1), nrow=3))
eigen(A, symmetric = TRUE)$values

Matrix::chol(A, pivot=FALSE) # 2nd argument requests pivoting

(Achol <- base::chol(A, pivot=TRUE)) # turn off checking pd
attr(Achol, "pivot")
attr(Achol, "rank")

Matrix::rankMatrix(A) # determine rank from SVD, which is more numerically stable

Achol[, attr(Achol, "pivot")]

(P <- as(attr(Achol, "pivot"), "pMatrix"))   # permutation matrix, actually P' in our notation

# P A P' = L L'
Matrix::norm(t(P) %*% A %*% P - t(Achol) %*% Achol, type="F")  # Frobenius norm

# word-by-word transcription of mathematical expression
# Matrix::determinant computes the logarithm of the determinant
logpdf_mvn_1 <- function(y, Sigma) {
    n <- length(y)
    - (n / 2) * log(2 * pi) - (1 / 2) * Matrix::determinant(Sigma)$modulus - (1 / 2) * t(y) %*% solve(Sigma) %*% y
}

# you learnt why you should avoid inversion
logpdf_mvn_2 <- function(y, Sigma) {
    n <- length(y)
    Sigmachol <- Matrix::chol(Matrix::symmpart(Sigma))
    - (n / 2) * log(2 * pi) - Matrix::determinant(Sigmachol)$modulus - (1 / 2) * sum(abs(solve(t(Sigmachol), y))^2)
}

# this is for the efficiency-concerned  
logpdf_mvn_3 <- function(y, Sigma) {
    n <- length(y)
    Sigmachol <- Matrix::chol(Matrix::symmpart(Sigma))
    - 0.5 * n * log(2 * pi) - sum(log(diag(Sigmachol))) - 0.5 * sum(abs(solve(t(Sigmachol), y))^2)
}

set.seed(123) # seed

n <- 1000
# a pd matrix
Sigma_half <- Matrix(rnorm(n * (2 * n)), nrow=n)
Sigma <- Sigma_half %*% t(Sigma_half)
Sigma <- Matrix::symmpart(Sigma)
Sigmachol <- Matrix::chol(Sigma)

y <- Sigmachol %*% rnorm(n)  # one random sample from N(0, Sigma). Homework

# at least they give same answer
(logpdf_mvn_1(y, Sigma))
(logpdf_mvn_2(y, Sigma))
(logpdf_mvn_3(y, Sigma))

res <- microbenchmark::microbenchmark(logpdf_mvn_1(y, Sigma), logpdf_mvn_2(y, Sigma), logpdf_mvn_3(y, Sigma))
print(res)