No pivoting for CHOLMOD factors

[jkbest2]

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.

[ElOceanografo]

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

[mschauer]

This bit me. As far as I can see, PDSparseMat seems to be non-functional and should not be offered to the ecosystem.