-
Notifications
You must be signed in to change notification settings - Fork 421
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
Comments
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 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? |
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.The text was updated successfully, but these errors were encountered: