Julia中Distributions.jl：如何用Cholesky矩阵高效定义MvNormal？

问题 #

我正在编写代码，主要通过操作协方差矩阵的Cholesky因子来变换多元分布。现有代码如下：

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

常规做法是：

c1 = cholesky(cov(g1)).L
# 对c1进行操作，示例中直接赋值给c2
c2 = c1
s2 = c2*c2'
g2 = MvNormal([1,2], s2)

但这种方式需要反复计算协方差矩阵的Cholesky因子再还原协方差矩阵，出于速度和数值稳定性考虑，我希望始终使用Cholesky因子（很少需要协方差矩阵）。即使简单示例也已出现数值误差：

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

我曾看到GitHub issue的相关回答提到可仅用Cholesky因子定义MvNormal，但因对Julia不够熟悉且方法已更新，无法实现。请问：若已有协方差矩阵的Cholesky因子，如何高效定义MvNormal，既能高效获取该因子，又能使用MvNormal的所有方法？注：c1到c2的变换是无迹变换平方根形式的实现，我希望自行实现以便调整。

解决方案 #

在Distributions.jl中，你可以直接基于Cholesky因子构造MvNormal，无需先还原协方差矩阵，既能避免数值误差，又能保留对Cholesky因子的直接访问。

两种高效构造方式 #

方式1：传入完整Cholesky对象

先构造出标准的Cholesky实例，直接作为协方差参数传入MvNormal：

using Distributions
g1 = MvNormal([1,2], [2 1; 1 2])
# 获取完整的Cholesky分解对象，而非仅下三角矩阵L
c1 = cholesky(cov(g1))
# 自行操作c1.L得到目标下三角因子c2_L（示例直接赋值）
c2_L = c1.L
# 基于新的L矩阵构造Cholesky对象，指定分解类型为下三角(:L)
c2 = Cholesky(c2_L, :L, 0)
# 直接用Cholesky对象构造多元正态分布
g2 = MvNormal([1,2], c2)

方式2：直接传入下三角因子L

Distributions.jl支持直接传入下三角矩阵作为协方差的平方根，内部会自动封装为Cholesky结构：

using Distributions
g1 = MvNormal([1,2], [2 1; 1 2])
c1_L = cholesky(cov(g1)).L
# 自行操作得到目标下三角因子c2_L
c2_L = c1_L
# 直接用L矩阵构造MvNormal
g2 = MvNormal([1,2], c2_L)

效果验证 #

用上述方式构造的g2不会出现数值误差，且后续调用cholesky(g2)时会直接返回存储的Cholesky因子，无需重新计算：

# 查看协方差矩阵，无数值偏差
cov(g2)
# 验证直接获取的Cholesky因子与操作后的L一致
cholesky(g2).L == c2_L  # 返回true

这种方式全程保留Cholesky因子的操作链路，避免了协方差矩阵的还原与重复分解，既提升了运算速度，又保证了数值稳定性，同时能正常使用MvNormal的所有内置方法（如抽样、概率密度计算等）。

内容的提问来源于stack exchange，提问作者Sirplentifus