C++ 实现 Wishart 分布样本矩阵生成

4.2 数值稳定性设计

1. 假设输入 Sigma 为对称正定
2. Cholesky 分解只使用下三角
3. 使用 double 精度
4. 明确矩阵乘法顺序减少误差

4.3 数据结构选择

• 矩阵：std::vector<std::vector<double>>
• 随机数：std::mt19937

五、完整实现代码

/******************************************************
 * File: wishart.h
 ******************************************************/
#ifndef WISHART_H
#define WISHART_H
#include <vector>

/* 从 Wishart 分布生成样本矩阵 */
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(
        const std::vector<std::vector<double>>& A);
    std::vector<std::vector<double>> multiply(
        const std::vector<std::vector<double>>& A,
        const std::vector<std::vector<double>>& B);
    std::vector<std::vector<double>> transpose(const std::vector<std::vector<double>>& A);
};
#endif

/******************************************************
 * File: wishart.cpp
 ******************************************************/
#include "wishart.h"
#include <cmath>
#include <random>
#include <stdexcept>

/* 构造函数 */
WishartSampler::WishartSampler(int dim, int df) : p(dim), n(df) {
    if (n < p) throw std::invalid_argument("df must be >= dimension");
}

/* Cholesky 分解（下三角） */
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;
}

/* 矩阵乘法 */
std::vector<std::vector<double>> WishartSampler::multiply(
    const std::vector<std::vector<double>>& A,
    const std::vector<std::vector<double>>& B) {
    std::vector<std::vector<double>> C(p, std::vector<double>(p, 0.0));
    for (int i = 0; i < p; ++i)
        for (int j = 0; j < p; ++j)
            for (int k = 0; k < p; ++k)
                C[i][j] += A[i][k] * B[k][j];
    return C;
}

/* 转置 */
std::vector<std::vector<double>> WishartSampler::transpose(
    const std::vector<std::vector<double>>& A) {
    std::vector<std::vector<double>> T(p, std::vector<double>(p));
    for (int i = 0; i < p; ++i)
        for (int j = 0; j < p; ++j)
            T[j][i] = A[i][j];
    return T;
}

/* Wishart 采样 */
std::vector<std::vector<double>> WishartSampler::sample(
    const std::vector<std::vector<double>>& sigma) {
    std::mt19937 gen(std::random_device{}());
    std::normal_distribution<> normal(0.0, 1.0);

    /* Step 1: Cholesky 分解 */
    auto L = choleskyDecomposition(sigma);

    /* Step 2: 构造 Bartlett 矩阵 A */
    std::vector<std::vector<double>> A(p, std::vector<double>(p, 0.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);
    }

    /* Step 3: B = L * A */
    auto B = multiply(L, A);

    /* Step 4: W = B * B^T */
    auto BT = transpose(B);
    return multiply(B, BT);
}

/******************************************************
 * File: main.cpp
 ******************************************************/
#include <iostream>
#include "wishart.h"

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 sample matrix:\n";
    for (const auto& row : W) {
        for (double v : row)
            std::cout << v << " ";
        std::cout << "\n";
    }
    return 0;
}

六、代码详细解读（仅解读方法作用）

6.1 sample

• 核心 Wishart 采样流程
• 实现 Bartlett 分解
• 返回对称正定矩阵

6.2 choleskyDecomposition

• 将尺度矩阵分解为下三角矩阵
• 引入协方差结构

6.3 Bartlett 矩阵构造

• 对角线使用卡方分布
• 下三角使用标准正态分布

七、项目详细总结

通过本项目，我们：

1. 系统理解了 Wishart 分布的定义与意义
2. 掌握了 Bartlett 分解的工程实现
3. 学会了随机矩阵的数值生成方法
4. 为贝叶斯协方差建模打下基础

该实现：

• 完全独立
• 数值稳定
• 教学清晰
• 可直接用于科研

八、项目常见问题及解答

Q1：为什么要求 df ≥ p？

A：这是 Wishart 分布存在的必要条件。

Q2：是否保证输出正定？

A：是的，Bartlett 分解理论上保证正定。

Q3：如何生成逆 Wishart 分布？

A：对 Wishart 矩阵取逆即可。

九、扩展方向与性能优化

1. 实现逆 Wishart 分布
2. 使用 Kahan 求和提高数值稳定性
3. 使用 SIMD 加速矩阵乘法
4. 扩展到稀疏 Wishart
5. 集成到贝叶斯 GMM