Open ForceBru opened 2 years ago
I used brute force to test whether the matrix is positive definite with the Sylvester's criterion:
using LinearAlgebra
gp_covar = kernel_rbf_mat(X, X)
determinants = @views [
det(gp_covar[begin:i, begin:i])
for i in axes(gp_covar, 1)
]
This gives me many zeros and some negative determinants that are really close to zero:
100-element Vector{Float64}:
1.0
0.010151166063869232
2.0814591040367652e-6
1.2865923518518836e-11
3.1803252426530575e-18
3.909951852971926e-26
2.8544529806688935e-35
1.4394891501197602e-45
6.387850389715162e-57
-1.5246070532959875e-68
3.315279979574806e-80
2.6577421180681176e-93
-2.437546507932742e-104
⋮
0.0
-0.0
-0.0
-0.0
-0.0
0.0
-0.0
-0.0
0.0
0.0
-0.0
0.0
So yes, the covariance matrix is not positive definite. However, NumPy's multivariate_normal
works fine anyway. In fact, its documentation says (emphasis mine):
Covariance matrix of the distribution. It must be symmetric and positive-semidefinite for proper sampling.
Wikipedia also says in the "Equivalent definitions" section:
A random vector ... has a multivariate normal distribution if ...:
There is a k-vector mu and a symmetric, positive semidefinite kxk matrix sigma...
So, looks like the covariance matrix should be positive semidefinite.
I'd say my covariance matrix is in fact positive semidefinite: numbers like -1.52e-68
, -2.4375e-104
and -0.0
are essentially zeros.
https://github.com/JuliaStats/Distributions.jl/issues/1602#issuecomment-1209901969 explains how one can use positive semi-definite matrices (to some extent, i.e., e.g. for rand
).
TL;DR: given the same mean vector and covariance matrix, NumPy can sample from a multivariate normal with these parameters, but
Distributions.MultivariateNormal
says:PosDefException: matrix is not positive definite; Cholesky factorization failed
.Julia code
Julia error
Apparently, the covariance matrix
gp_covar
is not positive definite:Equivalent Python code
Problem
How to reproduce:
julia-1.8 gp_sample.jl
. The covariance matrix will be saved tocovariance_jl.csv
.python3 gp_sample.py
. The covariance matrix will be saved tocovariance_py.csv
.However, given the same parameters, Python can sample from the multivariate normal no problem, but Distributions.jl can't.
Also, NumPy apparently doesn't care that the covariance matrix
gp_covar
isn't positive definite (is it actually not positive definite?? I'm not sure anymore) and samples anyway, while Distributions.jl thinks that it's a problem.This code is supposed to sample from a Gaussian process (here's the tutorial I'm following: https://peterroelants.github.io/posts/gaussian-process-tutorial/#Sampling-from-prior), so it... should work?
Simpler MWE
Actually, the kernel stuff isn't really needed to reproduce the bug. The covariance matrix can simply be loaded from
covariance_py.csv
:Versions