在C++中实现考虑自相关的多元似然函数：求协方差矩阵逆与行列式

对于你的多元似然函数需求，优先推荐使用成熟的线性代数库来处理矩阵运算——手动实现高阶矩阵的行列式和逆不仅容易出错，还会面临数值稳定性问题。下面分两种场景给出实现方案：

方案一：使用Eigen线性代数库（推荐） #

Eigen是C++生态中最常用的开源线性代数库，仅需头文件即可使用，支持高效的矩阵分解、行列式计算和求逆操作，尤其适合处理对称正定的协方差矩阵（考虑时间自相关的协方差矩阵通常满足正定条件）。

完整代码示例 #

#include <Eigen/Dense>
#include <cmath>
#include <iostream>

int main() {
    // 维度d：对应时间步或变量的数量
    int d = 3;

    // 1. 构造协方差矩阵Σ（替换为你的自相关协方差矩阵）
    Eigen::MatrixXd sigma(d, d);
    sigma << 1.0, 0.2, 0.1,
             0.2, 1.0, 0.3,
             0.1, 0.3, 1.0;

    // 2. 检查矩阵是否正定（协方差矩阵理论上正定，这里做容错）
    Eigen::LLT<Eigen::MatrixXd> llt_decomp(sigma);
    if (llt_decomp.info() != Eigen::Success) {
        std::cerr << "错误：协方差矩阵非正定或奇异，无法计算逆！" << std::endl;
        return 1;
    }

    // 3. 计算行列式|Σ|：利用Cholesky分解的L矩阵对角线元素平方乘积
    double det_sigma = std::pow(llt_decomp.matrixL().diagonal().prod(), 2);

    // 4. 计算逆矩阵Σ^{-1}：通过Cholesky分解求解单位矩阵得到
    Eigen::MatrixXd sigma_inv = llt_decomp.solve(Eigen::MatrixXd::Identity(d, d));

    // 5. 代入似然函数公式计算
    Eigen::VectorXd x(d);    // 模拟数据向量
    Eigen::VectorXd mu(d);   // 观测数据向量
    x << 2.1, 3.0, 1.8;
    mu << 2.0, 3.1, 1.7;
    Eigen::VectorXd diff = x - mu;

    // 计算二次型项 (x-μ)^T * Σ^{-1} * (x-μ)
    double quadratic_form = diff.transpose() * sigma_inv * diff;

    // 计算似然函数值
    double likelihood = 1.0 / std::sqrt(std::pow(2 * M_PI, d) * det_sigma) * std::exp(-0.5 * quadratic_form);

    std::cout << "似然函数值：" << likelihood << std::endl;

    return 0;
}

关键说明 #

• Cholesky分解（LLT类）是处理对称正定矩阵的最优选择，比通用逆运算inverse()效率更高、数值更稳定；
• 必须检查矩阵是否正定/非奇异，避免求逆时出现错误；
• Eigen的MatrixXd和VectorXd支持动态维度，可直接适配你的实际数据规模。

方案二：手动实现（仅用于学习原理） #

如果想理解底层逻辑，可以手动实现小维度矩阵的行列式和逆计算，但不推荐用于实际项目（高阶矩阵会出现数值误差和效率问题）。

3x3矩阵的手动实现示例 #

#include <cmath>
#include <iostream>
#include <vector>
#include <stdexcept>

using Matrix = std::vector<std::vector<double>>;
using Vector = std::vector<double>;

// 计算3x3矩阵的行列式
double calc_det_3x3(const Matrix& mat) {
    return mat[0][0] * (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1])
         - mat[0][1] * (mat[1][0] * mat[2][2] - mat[1][2] * mat[2][0])
         + mat[0][2] * (mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0]);
}

// 计算3x3矩阵的逆
Matrix calc_inv_3x3(const Matrix& mat) {
    double det = calc_det_3x3(mat);
    if (std::abs(det) < 1e-9) {
        throw std::runtime_error("奇异矩阵无法求逆");
    }

    Matrix inv(3, std::vector<double>(3));
    inv[0][0] = (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]) / det;
    inv[0][1] = (mat[0][2] * mat[2][1] - mat[0][1] * mat[2][2]) / det;
    inv[0][2] = (mat[0][1] * mat[1][2] - mat[0][2] * mat[1][1]) / det;
    inv[1][0] = (mat[1][2] * mat[2][0] - mat[1][0] * mat[2][2]) / det;
    inv[1][1] = (mat[0][0] * mat[2][2] - mat[0][2] * mat[2][0]) / det;
    inv[1][2] = (mat[0][2] * mat[1][0] - mat[0][0] * mat[1][2]) / det;
    inv[2][0] = (mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0]) / det;
    inv[2][1] = (mat[0][1] * mat[2][0] - mat[0][0] * mat[2][1]) / det;
    inv[2][2] = (mat[0][0] * mat[1][1] - mat[0][1] * mat[1][0]) / det;

    return inv;
}

// 计算二次型项
double calc_quadratic_form(const Vector& diff, const Matrix& inv_mat) {
    double result = 0.0;
    int d = diff.size();
    for (int i = 0; i < d; ++i) {
        double row_sum = 0.0;
        for (int j = 0; j < d; ++j) {
            row_sum += diff[j] * inv_mat[i][j];
        }
        result += diff[i] * row_sum;
    }
    return result;
}

int main() {
    try {
        int d = 3;
        Matrix sigma = {{1.0, 0.2, 0.1}, {0.2, 1.0, 0.3}, {0.1, 0.3, 1.0}};
        double det_sigma = calc_det_3x3(sigma);
        Matrix sigma_inv = calc_inv_3x3(sigma);

        Vector x = {2.1, 3.0, 1.8};
        Vector mu = {2.0, 3.1, 1.7};
        Vector diff(d);
        for (int i = 0; i < d; ++i) {
            diff[i] = x[i] - mu[i];
        }

        double quadratic_form = calc_quadratic_form(diff, sigma_inv);
        double likelihood = 1.0 / std::sqrt(std::pow(2 * M_PI, d) * det_sigma) * std::exp(-0.5 * quadratic_form);

        std::cout << "似然函数值：" << likelihood << std::endl;
    } catch (const std::exception& e) {
        std::cerr << "错误：" << e.what() << std::endl;
        return 1;
    }

    return 0;
}

关键说明 #

• 手动实现仅适合小维度矩阵（如3x3），高阶矩阵需用LU分解等更复杂的方法；
• 必须检查行列式是否接近0，避免除以0的错误；
• 手动实现的数值稳定性远低于成熟库，实际项目中不建议使用。

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