The covariance matrix of my multivariate normal is not full rank and so regular cholesky throws a positive definite error.
However, it is possible to do a pivoted cholesky decomposition. Quoting from the netlib manual of pstrf : http://www.netlib.org/lapack/explore-html/dd/dad/dpstrf_8f.html The factorization has the form P**T * A * P = U**T * U , if UPLO = 'U', P**T * A * P = L * L**T, if UPLO = 'L', where U is an upper triangular matrix and L is lower triangular, and P is stored as vector PIV. I want to generate a MvNormal draw by multiplying P*L*z where z is a vector of standard normals. >From the lapack doc it seems I should be able to recover A from P*L*L'*P'. I have a minimial reproducible example where I try to recover A, but some of the elements of the difference matrix are larger than i would expect or hope: srand(1) a=rand(3,10) A=a'*a rank(A) F=cholfact(A,:L,Val{true}) L=F[:L] P=F[:P] P*L*L'*P'-A julia> P*L*L'*P'-A 10x10 Array{Float64,2}: 1.07228 0.0 0.0 1.27149 1.11022e-16 0.0 -5.55112e-17 1.84107 2.01638 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -2.22045e-16 0.0 0.0 0.0 2.22045e-16 0.0 0.0 1.27149 0.0 -2.22045e-16 3.11518 -3.33067e-16 0.0 0.0 2.82245 2.48629 0.0 1.11022e-16 0.0 0.0 -3.33067e-16 5.55112e-16 0.0 0.0 0.0 2.22045e-16 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -5.55112e-17 0.0 0.0 0.0 -5.55112e-17 0.0 0.0 0.0 0.0 -5.55112e-17 0.0 0.0 0.0 0.0 1.84107 0.0 2.22045e-16 2.82245 0.0 0.0 0.0 4.15206 3.59957 0.0 2.01638 0.0 0.0 2.48629 2.22045e-16 0.0 0.0 3.59957 4.09992 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 notice some elements are large. Why is it not identically zero?
