A = [4 12 -16; 12 37 -43; -16 -43 98]

# Cholesky without pivoting
Achol = cholfact(A)

Achol[:L]

Achol[:U]

b = [1.0; 2.0; 3.0]
A \ b # this does LU; wasteful!; 2/3 n^3 + 2n^2

@which A \ b

Achol \ b # two triangular solves; only 2n^2 flops

@which Achol \ b

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

@which det(A)

det(Achol) # cheap

@which det(Achol)

inv(A) # this does LU!

@which inv(A)

inv(Achol) 

@which inv(Achol)

srand(280) # seed
A = randn(5, 3)
A = A * A'

Achol = cholfact(A)

Achol = cholfact(A, :L, Val{true}) # 3rd argument request partial pivoting

rank(Achol) # determine rank from Cholesky factor

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

Achol[:L]

Achol[:U]

Achol[:p]

# P A P' = L U
vecnorm(Achol[:P] * A * Achol[:P]' - Achol[:L] * Achol[:U])

# this is a person w/o any numerical analsyis training
function logpdf_mvn_1(y::Vector, Σ::Matrix)
    n = length(y)
    - (n//2) * log(2π) - (1//2) * logdet(Σ) - (1//2) * y' * inv(Σ) * y
end

# this is an efficiency-savvy person 
function logpdf_mvn_2(y::Vector, Σ::Matrix)
    n = length(y)
    Σchol = cholfact(Symmetric(Σ))
    - (n//2) * log(2π) - (1//2) * logdet(Σchol) - (1//2) * sum(abs2, Σchol[:L] \ y)
end

# better memory efficiency
function logpdf_mvn_3(y::Vector, Σ::Matrix)
    n = length(y)
    Σchol = cholfact(Symmetric(Σ))
    - (n//2) * log(2π) - (1//2) * logdet(Σchol) - (1//2) * dot(y, Σchol \ y)
end

using BenchmarkTools, Distributions

srand(280) # seed

n = 1000
Σ = randn(n, n); Σ = Σ' * Σ + I # covariance matrix
y = rand(MvNormal(Σ)) # one randdom sample from N(0, Σ)

# at least they give same answer
@show logpdf_mvn_1(y, Σ)
@show logpdf_mvn_2(y, Σ)
@show logpdf_mvn_3(y, Σ);

@benchmark logpdf_mvn_1(y, Σ)

@benchmark logpdf_mvn_2(y, Σ)

@benchmark logpdf_mvn_3(y, Σ)

versioninfo()