invenia / PDMatsExtras.jl

Extra Positive (Semi-)Definite Matricies
MIT License
8 stars 6 forks source link

Fix pivoted order in whiten & unwhiten #33

Closed vandenman closed 1 year ago

vandenman commented 1 year ago

This bug popped up in Slack when rand(Distributions.MvNormal(PSDMat(...))) seemingly returned permuted samples. This PR is a proposal to fix that.

on master

using PDMats, PDMatsExtras

mat = [1.; 0 ;; 0 ; 2]
pd  = PDMat(mat)
psd = PSDMat(mat)
psd.chol.p # [2, 1]

x = ones(2)
whiten!(similar(x), pd,  x)
whiten!(similar(x), psd, x)

x = ones(2)
unwhiten!(similar(x), pd,  x)
unwhiten!(similar(x), psd, x)

x = ones(2, 3)
whiten!(similar(x), pd,  x) 
whiten!(similar(x), psd, x)

x = ones(2, 3)
unwhiten!(similar(x), pd,  x) 
unwhiten!(similar(x), psd, x)

For PDMat this returns e.g.,

2-element Vector{Float64}:
 1.0
 0.7071067811865475

but for PSDMat this returns

2-element Vector{Float64}:
 0.7071067811865475
 1.0

because the pivot in psd.chol.p is ignored.

This PR fixes this and adds a test.

Comments, feedback, and other suggestions are welcome!