Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Closed
jkbest2 opened this issue Oct 16, 2020 · 1 comment · Fixed by #1218
Closed

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

jkbest2 opened this issue Oct 16, 2020 · 1 comment · Fixed by #1218

Comments

@jkbest2
Copy link
Contributor

jkbest2 commented Oct 16, 2020

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
Copy link
Contributor

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"))

image

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"))

image

@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?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants