使用Cholesky分解生成多元正态矩阵的精度问题

我之前在手动实现多元正态采样的时候也碰到过一模一样的问题——当均值向量U的量级远小于X@L输出的波动幅度时，直接相加会因为浮点数精度限制导致U的有效位数被覆盖，最终模拟结果的均值和预期偏差明显。下面给你几个实用的解决思路：

• 先生成零均值样本，再单独平移均值
  这是最直接的解决方案：先基于Cholesky分解生成零均值的多元正态样本Y0 = X@L，然后逐维度单独加上对应的U分量，而不是把U和矩阵运算放在一起。这样能避免浮点数相加时，小量级的U被大量级的X@L结果“吃掉”。举个Python示例：

  import numpy as np
  n, m = 1000, 5
  X = np.random.normal(size=(n, m))  # n×m标准正态矩阵
  Σ = np.random.rand(m, m)
  Σ = Σ @ Σ.T  # 构造正定协方差矩阵
  L = np.linalg.cholesky(Σ)
  U = np.array([1e-8, 2e-8, 3e-8, 4e-8, 5e-8])  # 接近0的均值向量
  
  # 正确的做法：先生成零均值样本，再逐列加均值
  Y0 = X @ L
  Y = Y0 + U.reshape(1, m)  # 确保U被广播到每一行
  
  # 验证均值
  print(np.mean(Y, axis=0))  # 现在应该更接近U的真实值
  

• 调整变量缩放比例
  如果U接近0是因为协方差矩阵的尺度太大，你可以先对协方差矩阵做缩放，生成样本后再反缩放并加上均值：

  1. 选一个合适的缩放因子s（比如让缩放后的协方差矩阵的对角线元素量级接近1）
  2. 计算缩放后的协方差矩阵Σ_scaled = Σ / s²，对应的Cholesky因子L_scaled = L / s
  3. 生成零均值样本Y0_scaled = X @ L_scaled
  4. 还原尺度并加均值：Y = Y0_scaled * s + U
     这样X@L_scaled的量级会变小，U的相对量级就提升了，精度损失会大幅减少。

• 切换到更高精度的数值类型
  如果你的编程语言支持，可以把所有计算切换到更高精度的浮点数类型（比如Python的np.float128、R的longdouble）。更高的精度能保留更多有效位数，即使U的量级很小，也不会被X@L的结果覆盖。示例：

  # 切换到float128精度
  X = np.random.normal(size=(n, m)).astype(np.float128)
  L = np.linalg.cholesky(Σ).astype(np.float128)
  U = U.astype(np.float128)
  Y = X @ L + U
  

  注意：更高精度的计算会有一定性能开销，对于大规模数据集需要权衡速度和精度。

• 直接使用内置的多元正态生成函数
  其实大多数科学计算库（比如NumPy的np.random.multivariate_normal、R的MASS::mvrnorm）已经内置了稳定的多元正态采样实现，它们内部已经处理了小均值的精度问题。如果不是必须手动实现Cholesky分解，直接用这些内置函数会更省心，也能避免手动实现带来的数值问题。示例：

  Y = np.random.multivariate_normal(mean=U, cov=Σ, size=n)
  

总结一下：如果不需要手动实现Cholesky逻辑，优先用内置函数；如果必须自己实现，先生成零均值样本再单独加均值是最优选择；调整缩放比例和使用高精度类型是备选方案。

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