在 Julia 的 Distributions.jl 包中，如何使用 Cholesky 矩阵定义 MvNormal 分布？

Question

我正在编写一些代码，主要通过对协方差矩阵的胆汁因子进行操作来转换多变量分布。假设我有

using Distributions
g1 = MvNormal([1,2], [2 1; 1 2])

天真的做法是

c1 = cholesky(cov(g1)).L
# work on matrix c1 in some way, for this example I'll just directly assign it to c2 
c2 = c1
s2 = c2*c2'
g2 = MvNormal([1,2], s2)

但这意味着重复计算协方差矩阵的乔列斯基因子，然后获得相应的协方差矩阵，而出于速度和数值稳定性的原因，我更愿意始终使用乔列斯基因子，因为我很少需要协方差矩阵。请注意，即使使用这个非常简单的示例，也可以看到一些错误：

g2 = MvNormal([1,2], s2)
FullNormal(
dim: 2
μ: [1.0, 2.0]
Σ: [2.0000000000000004 1.0; 1.0 1.9999999999999996]
)

Julia 执行此操作的正确方法是什么？例如，我见过 GitHub 问题的答案，它表明可以仅指定 MvNormal 的 Cholesky 因子，但它假设我对 Julia 有一定的熟悉度，而我并不熟悉，所以这对我来说还不够弄清楚如何实际做到这一点。而且它已经很老了，周围的讨论表明它的工作方式同时已经改变了。

因此，总而言之，如果我有协方差矩阵的胆汁因子，我如何有效地定义 MvNormal，以便我可以有效地获取该因子，并且仍然可以访问为 MvNormal 定义的方法？

注意：我的应用程序中

c1

到

c2

的变换是无迹变换的平方根形式的实现。我意识到可能有库可以做到这一点，但我更愿意有自己的实现，因为我计划修改它。

Answer 1

使用 PDMats 包我想你可以获得你需要的功能：

julia> using Distributions, LinearAlgebra, PDMats

julia> R = rand(5,5);

julia> M = R*transpose(R);

julia> C = cholesky(M)
Cholesky{Float64, Matrix{Float64}}
U factor:
5×5 UpperTriangular{Float64, Matrix{Float64}}:
 1.44058  1.09028   0.971589  1.09389   1.36865
  ⋅       0.939203  0.783269  0.721269  0.0734465
  ⋅        ⋅        0.576735  0.286545  0.614992
  ⋅        ⋅         ⋅        0.892629  0.662157
  ⋅        ⋅         ⋅         ⋅        0.256624

julia> N = PDMat(C)
5×5 PDMat{Float64, Matrix{Float64}}:
 2.07526  1.57064  1.39965  1.57583  1.97165
 1.57064  2.07082  1.79496  1.87007  1.5612
 1.39965  1.79496  1.89012  1.79302  1.74198
 1.57583  1.87007  1.79302  2.59572  2.31741
 1.97165  1.5612   1.74198  2.31741  2.76113

julia> N ≈ M
true

julia> mvn = MvNormal(zeros(5), N);

julia> rand(mvn, 2)
5×2 Matrix{Float64}:
 -3.64747   -1.30948
 -1.84731   -1.00835
 -1.68674   -1.01449
 -0.818989   2.4971
 -2.87103    1.37066

所以现在，

mvn

应该是 MvNormal，其中

C

作为 Cholesky 因子。任何新的 Cholesky 因子都可以通过从 Cholesky 因子对象构造

PDMat

（可以从下三角矩阵构造）来以这种方式生成 MvNormal。

希望开销小于每次生成新的 MvNormal 协方差矩阵。