PDMats doesn't account for pivoting in Cholesky of sparse matrices

[jkbest2]

See JuliaStats/PDMats.jl#120.

This results in incorrect realizations from rand(::MvNormalCanon) when the provided precision matrix is factored by CHOLMOD. It may affect other methods as well.

[ElOceanografo]

MWE of this with a Gaussian Markov random field on a grid:

using Distributions, PDMats
using SparseArrays, LinearAlgebra
using Plots 

r = 0.25
n = 50
J = spdiagm(-n => fill(-r, n*(n-1)),
            -1 => fill(-r, n^2 - 1), 
            0 => fill(1, n^2),
            1 => fill(-r, n^2 - 1),
            n => fill(-r, n*(n-1)))
dsparse = MvNormalCanon(PDSparseMat(J))
ddense = MvNormalCanon(Matrix(J))

x1 = rand(dsparse)
x2 = rand(ddense)
plot(heatmap(reshape(x1, n, n), title="Sparse"), heatmap(reshape(x2, n, n), title="Dense"))

Fortunately, the pivoting doesn't affect likelihood calculations:

logpdf(dsparse, x1) ≈ logpdf(ddense, x1) # true

As I wrote in #1201, this is an easy code fix, but it requires thinking about how to rewrite some of the tests (which are probably broken currently, see this comment on #855). The Cholesky decomposition is not unique, and the implementations in LinearAlgebra/LAPACK and SuiteSparse/CHOLMOD appear to return different factorizations. As a result, the samples obtained from a dense and sparse MvNormalCanon will be different, even if their precision matrices and the random seeds are identical. (The samples do have the same distribution, however.)

z = randn(n^2)
x3 = dsparse.J.chol.PtL' \ z
x4 = ddense.J.chol.U \ z
plot(heatmap(reshape(x3, n, n), title="Sparse"), heatmap(reshape(x4, n, n), title="Dense"))

@matbesancon @andreasnoack how do you feel about this? How much of a problem is it if the particular samples drawn from a distribution depend on the types of its parameters?