1. 项目概述
在多元统计分析和机器学习领域，Wishart分布是一个极其重要的矩阵值概率分布。作为一名长期从事统计计算和机器学习算法开发的工程师，我经常需要在项目中实现各种概率分布的采样功能。今天我要分享的是如何在C++中从零实现Wishart分布的采样器，不依赖任何第三方数学库。
Wishart分布可以看作是多元正态分布样本协方差矩阵的分布。简单来说，如果我们有n个来自p维多元正态分布的独立样本，那么这些样本的协方差矩阵就服从Wishart分布。这个分布在贝叶斯统计、金融风险建模和高斯过程等领域都有广泛应用。
2. Wishart分布理论基础
2.1 数学定义
Wishart分布W(n, Σ)的概率密度函数为：
f(W) = [|W|^{(n-p-1)/2} exp(-tr(Σ^{-1}W)/2)] / [2^{np/2} |Σ|^{n/2} Γ_p(n/2)]
其中：
- n是自由度（必须≥p）
- Σ是p×p的尺度矩阵（对称正定）
- Γ_p是多元Gamma函数
- tr表示矩阵迹
- |·|表示矩阵行列式
2.2 采样原理
从Wishart分布采样的核心方法是Bartlett分解法，其步骤如下：
- 对尺度矩阵Σ进行Cholesky分解：Σ = LLᵀ
- 生成Bartlett矩阵A：
- 对角线元素：A_{ii} ~ √χ²(n-i+1)
- 下三角元素：A_{ij} ~ N(0,1) (i>j)
- 上三角元素：A_{ij} = 0 (i<j)
- 计算B = LA
- 返回W = BBᵀ ~ W(n, Σ)
这种方法的优势在于数值稳定性好，且不需要依赖复杂的矩阵运算库。
3. C++实现详解
3.1 类结构设计
我们设计一个WishartSampler类，主要包含以下成员：
cpp复制class WishartSampler {
public:
WishartSampler(int dim, int df); // 构造函数
std::vector<std::vector<double>> sample(
const std::vector<std::vector<double>>& sigma); // 采样函数
private:
int p; // 维度
int n; // 自由度
// 辅助函数
std::vector<std::vector<double>> choleskyDecomposition(...);
std::vector<std::vector<double>> multiply(...);
std::vector<std::vector<double>> transpose(...);
};
3.2 Cholesky分解实现
Cholesky分解是算法的关键步骤之一，我们将Σ分解为下三角矩阵L：
cpp复制std::vector<std::vector<double>>
WishartSampler::choleskyDecomposition(
const std::vector<std::vector<double>>& A)
{
std::vector<std::vector<double>> L(p,
std::vector<double>(p, 0.0));
for (int i = 0; i < p; ++i) {
for (int j = 0; j <= i; ++j) {
double sum = 0.0;
for (int k = 0; k < j; ++k)
sum += L[i][k] * L[j][k];
if (i == j)
L[i][j] = std::sqrt(A[i][i] - sum);
else
L[i][j] = (A[i][j] - sum) / L[j][j];
}
}
return L;
}
注意：输入矩阵A必须是对称正定的，否则分解会失败。在实际应用中应该添加校验逻辑。
3.3 Bartlett矩阵构造
Bartlett矩阵的构造需要生成卡方分布和正态分布的随机数：
cpp复制// 构造Bartlett矩阵A
std::vector<std::vector<double>> A(p,
std::vector<double>(p, 0.0));
std::mt19937 gen(std::random_device{}());
std::normal_distribution<> normal(0.0, 1.0);
for (int i = 0; i < p; ++i) {
// 对角线元素：卡方分布
double chi2 = 0.0;
for (int k = 0; k < n - i; ++k) {
double z = normal(gen);
chi2 += z * z;
}
A[i][i] = std::sqrt(chi2);
// 下三角元素：标准正态分布
for (int j = 0; j < i; ++j)
A[i][j] = normal(gen);
}
3.4 完整采样流程
将各个步骤组合起来完成采样：
cpp复制std::vector<std::vector<double>>
WishartSampler::sample(
const std::vector<std::vector<double>>& sigma)
{
// Step 1: Cholesky分解
auto L = choleskyDecomposition(sigma);
// Step 2: 构造Bartlett矩阵A
// ... (见3.3节代码)
// Step 3: B = L * A
auto B = multiply(L, A);
// Step 4: W = B * B^T
auto BT = transpose(B);
return multiply(B, BT);
}
4. 数值稳定性与优化
4.1 数值稳定性考虑
- 输入验证：确保n ≥ p，Σ对称正定
- 双精度浮点：使用double而非float减少舍入误差
- 分解稳定性：Cholesky分解中添加小量正则化(Σ + εI)可增强稳定性
- 乘法顺序：矩阵乘法时合理安排计算顺序减少误差累积
4.2 性能优化建议
- 内存布局：使用一维数组存储矩阵，提高缓存利用率
- 并行计算：矩阵乘法可并行化
- SIMD指令：使用AVX等指令集加速计算
- 随机数生成：使用更高效的随机数生成器
5. 使用示例
下面是一个完整的使用示例：
cpp复制int main() {
int p = 3; // 维度
int df = 6; // 自由度
// 尺度矩阵（必须对称正定）
std::vector<std::vector<double>> sigma = {
{1.0, 0.5, 0.2},
{0.5, 1.0, 0.3},
{0.2, 0.3, 1.0}
};
WishartSampler sampler(p, df);
auto W = sampler.sample(sigma);
std::cout << "Wishart样本矩阵:\n";
for (const auto& row : W) {
for (double v : row)
std::cout << v << " ";
std::cout << "\n";
}
return 0;
}
6. 常见问题与解决方案
6.1 采样结果不是正定矩阵
可能原因：
- 自由度n太小（应满足n ≥ p）
- 尺度矩阵Σ不是正定的
- 数值误差累积
解决方案：
- 检查输入参数
- 对Σ进行特征值分解验证
- 添加小的正则化项
6.2 性能瓶颈
主要瓶颈通常在矩阵乘法部分。优化建议：
- 使用更高效的矩阵乘法实现
- 考虑使用Strassen算法
- 对于大矩阵，使用分块计算
6.3 扩展到高维情况
当维度p很大时：
- 考虑稀疏矩阵表示
- 使用近似方法
- 利用矩阵的特殊结构（如对角占优）
7. 扩展应用
7.1 逆Wishart分布
逆Wishart分布在贝叶斯统计中常用作协方差矩阵的先验。可以通过对Wishart样本取逆得到：
cpp复制auto inv_W = inverse(W); // 需要实现矩阵求逆
7.2 贝叶斯应用
在贝叶斯线性回归中，可以使用Wishart分布作为协方差矩阵的先验：
code复制p(Σ) ~ W(ν₀, S₀)
p(β|Σ) ~ N(μ₀, Σ⊗(X'X)⁻¹)
7.3 金融风险建模
在投资组合优化中，Wishart分布可用于建模资产收益的协方差矩阵不确定性。
8. 工程实践建议
在实际项目中应用时，我有以下几点建议：
- 单元测试：为关键函数（如Cholesky分解）编写详尽的测试用例
- 性能分析：使用profiler识别热点代码
- 接口设计：考虑支持多种矩阵类型（Eigen, Armadillo等）
- 异常处理：完善错误处理机制
- 文档注释：为每个函数添加详细的API文档
我在实际项目中遇到过的一个典型问题是数值稳定性。当矩阵条件数很大时，简单的Cholesky分解可能会失败。解决方案是添加一个小的正则化参数：
cpp复制// 正则化Cholesky分解
for (int i = 0; i < p; ++i) {
A[i][i] += 1e-10; // 添加小量
}
另一个实用技巧是预先分配所有内存，避免在循环中频繁申请释放内存，这可以显著提高性能。