Sampling from a dirac delta

[cstjean]

I occasionally like having Gaussians with covariance=0 (this is natural, given the formulation in MiniKalman.jl). The math for Kalman filtering works fine, but when sampling I have an issue, which is exacerbated by https://github.com/JuliaArrays/StaticArrays.jl/pull/533:

julia> rand(Gaussian([1.0,2], SDiagonal(0.0,0.0)))   # before StaticArrays#533
2-element Array{Float64,1}:
 1.0
 2.0

julia> rand(Gaussian([1.0,2], Diagonal(SVector(0.0,0.0))))
ERROR: PosDefException: matrix is not positive definite; Cholesky factorization failed.
Stacktrace:
 [1] #cholesky!#127(::Bool, ::Function, ::Diagonal{Float64,MArray{Tuple{2},Float64,1,2}}, ::Val{false}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/diagonal.jl:524
 [2] #cholesky#128 at ./none:0 [inlined]
 [3] cholesky at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/diagonal.jl:533 [inlined] (repeats 2 times)
 [4] rand(::Random.MersenneTwister, ::Gaussian{Array{Float64,1},Diagonal{Float64,SArray{Tuple{2},Float64,1,2}}}) at /home/cst-jean/.julia/packages/GaussianDistributions/k5dmO/src/GaussianDistributions.jl:82
 [5] rand(::Gaussian{Array{Float64,1},Diagonal{Float64,SArray{Tuple{2},Float64,1,2}}}) at /home/cst-jean/.julia/packages/GaussianDistributions/k5dmO/src/GaussianDistributions.jl:80

cholesky of an all-0 Diagonal matrix fails. I don't know where/if it should be special-cased. Any thought? According to Wikipedia, there can be a cholesky decomposition of a positive semi-definite matrix, but that case doesn't seem to be handled by LinearAlgebra.

[mschauer]

Related to e.g. #27249. There is check=false we can try, but maybe this gives garbage instead of positive semidefinite factors. In the end maybe we have to diverge from Base and provide our own chol wrapper here I think.

[mschauer]

You can of course use PSD and provide PSD(sqrt(Diagonal(v))) but further on you might need to define some methods as +(a::PSD, b::PSD) = PSD(sqrt((a.σ*a.σ') + (b.σ*b.σ')) etc. I think this would have the side effect of implementing the square root Kalman filter...