PDMats.jl
PDMats.jl copied to clipboard
No pivoting for CHOLMOD factors
whiten!
and unwhiten
use chol_lower
to extract the lower triangular component of Cholesky decompositions. For SuiteSparse.CHOLMOD.Factor{T}
, this extracts the lower triangular factor, but does not account for pivoting. This results in the same problem seen on this Discourse thread.
Using the sparse precision matrix attached to that thread,
using DataFrames, CSV
q_df = CSV.read("precision_matrix.csv")
point_df = CSV.read("mesh_points.csv")
Q = sparse(q_df.i, q_df.j, q_df.q)
Generate a realization,
using Random, SparseArrays, LinearAlgebra
using Distributions
using PDMats
Random.seed!(2)
z = randn(size(Q, 1))
pdQ = PDSparseMat(cholesky(Q))
xldiv = pdQ.chol.UP \ z
xwh = whiten(pdQ, z)
These two realizations should be the same, but they are not.
julia> maximum(xldiv - xwh)
45.00349329625478
This is also clear when they are plotted.
using Plots
pldiv = surface(point_df.x, point_df.y, xqs, camera=(45, 80), title="ldiv")
pwh = surface(point_df.x, point_df.y, xwh, camera=(45, 80), title="whiten")
plot(pldiv, pwh, layout=(1,2), size=(1200, 400), colorbar=false)
The ldiv
version results in a random field with the expected structure, while the whiten
version does not.
Downstream, this means that at least rand(::Distributions.MvNormalCanon)
is giving incorrect results. I haven't looked closely at other methods in Distributions.jl yet.
This issue may be worth fixing in PDMats, but note that the unwhiten_winv!
method for rand(::Distributions.MvNormalCanon)
currently extracts the Cholesky factor directly, without calling chol_lower
:
https://github.com/JuliaStats/Distributions.jl/blob/50a712acb24cedda4291e26bdf52cd7f2096da89/src/multivariate/mvnormalcanon.jl#L179
This bit me. As far as I can see, PDSparseMat
seems to be non-functional and should not be offered to the ecosystem.