如何在C语言中使用GSL的Cholesky分解函数获取下三角矩阵

首先纠正一个关键错误：你写的gsl_linalg_cholesky_decomp1不是GSL的标准函数，正确的Cholesky分解函数是gsl_linalg_cholesky_decomp。这个函数的核心特性是直接修改传入的矩阵，分解完成后，原矩阵的**下三角区域（含对角线）**就是你需要的Cholesky因子L（满足原矩阵A = L * L^T），上三角区域的值会保留原矩阵的内容，可以忽略。

以下是完整的实现代码：

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <stdio.h>

#define rows 6
#define cols 6

int main() {
    double cov[rows*cols] = {
        107.3461, 12.0710, -48.3746, 174.7796, 21.0202, -80.6075,
        12.0710, 8.0304, -5.9610, 20.2434, 2.2427, -9.3120,
        -48.3746, -5.9610, 25.2222, -78.6277, -9.4400, 36.1789,
        174.7796, 20.2434, -78.6277, 291.3491, 35.0176, -134.3626,
        21.0202, 2.2427, -9.4400, 35.0176, 4.2144, -16.1499,
        -80.6075, -9.3129, 36.1789, -134.3626, -16.1499, 61.9666
    };

    // 创建矩阵视图
    gsl_matrix_view m = gsl_matrix_view_array(cov, rows, cols);
    gsl_matrix *mat = &m.matrix;

    // 执行Cholesky分解
    int status = gsl_linalg_cholesky_decomp(mat);
    if (status != GSL_SUCCESS) {
        fprintf(stderr, "Cholesky分解失败，矩阵可能不是对称正定的\n");
        return 1;
    }

    // 提取并打印下三角矩阵L
    printf("Cholesky分解得到的下三角矩阵L:\n");
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < cols; j++) {
            if (i >= j) {
                // 下三角区域（含对角线）是L的元素
                printf("%10.4f ", gsl_matrix_get(mat, i, j));
            } else {
                // 上三角区域可忽略，这里打印0占位
                printf("%10.4f ", 0.0);
            }
        }
        printf("\n");
    }

    return 0;
}

关键说明： #

• gsl_linalg_cholesky_decomp会直接修改传入的矩阵，不需要额外分配内存存储结果，节省资源。
• 只有当原矩阵是对称正定矩阵时，分解才会成功，你的协方差矩阵满足这个条件，所以可以正常执行。
• 如果需要单独存储L矩阵，可以创建一个新的gsl_matrix，然后只复制原矩阵的下三角元素过去。

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