如何优化Cholesky分解求解的C语言实现以提升运行速度

问题描述 #

我希望使用Cholesky分解求解矩阵方程Ax=B，其中对称稠密矩阵A可分解为L和L^{T，转而求解L(L}Tx)=B：先通过前向替换求解Ly=B，再通过后向替换求解L^Tx=y。我用C语言实现了该逻辑，但与使用numpy.linalg.cho_solve和linalg.cho_factor的Python实现相比，C代码运行速度明显更慢。请问如何优化我的C代码以提升运行速度？

我的C代码实现 #

Cholesky分解函数 #

double* cholesky_factor(double* matrix, size_t n_rows, size_t n_cols) {
    
    // Timer
    clock_t tik, tok;
    double time_taken;
    tik = clock();

    double* L = (double*)calloc(n_rows * n_cols, sizeof(double));
    if (L == NULL) {
        fprintf(stderr, "Error: calloc failed to allocate memory\n");
        exit(EXIT_FAILURE);
    }
    
    printf("start");

    for(size_t i = 0; i < n_rows; i++){
        for (size_t j = 0; j <= i; j++){
            double sum = 0;
            for (size_t k = 0; k < j; k++)
            {
                sum += L[i * n_cols + k] * L[j * n_cols + k];
            }
            
            if (i == j) {
                L[i * n_cols + i] = sqrt(matrix[i * n_cols + i] - sum);
            } else {
                L[i * n_cols + j] = (1.0 / L[j * n_cols + j]) * (matrix[i * n_cols + j] - sum);
            }
        }
    }
        
    tok = clock();
    time_taken = (double)(tok - tik) / CLOCKS_PER_SEC;
    printf("C timer, Time taken for cholesky_factor: %f seconds\n", time_taken);

    return L;
}

基于Cholesky因子的求解函数 #

void transpose(int n_rows, int n_cols, const double* matrix, double* tran) {

    for (int i = 0; i < n_rows; i++) {
        for (int j = 0; j < n_cols; j++) {
            tran[j * n_rows + i] = matrix[i * n_cols + j];
        }
    }
}

double* forward_substitusion(const double* L, const double* B, size_t n_rows,size_t n_cols){
    double* Y = (double*)calloc(n_rows * n_cols, sizeof(double));
    
    for (size_t i = 0; i < n_rows; i++) {
        for (size_t j = 0; j < n_cols; j++) {
            double sum = 0;
            for (size_t k = 0; k < n_rows; k++) {
                sum += L[i*n_rows + k] * Y[k*n_rows + j];
            }
            Y[i*n_rows + j] = (B[i*n_rows + j] - sum) / L[i*n_rows + i];
        }
    }

    return Y;
}

double* backwards_substitution(const double* U, const double* y, size_t n_rows, size_t n_cols) {
    double* X = (double*)calloc(n_rows * n_cols, sizeof(double));

    for (int i = n_rows - 1; i >= 0; i--) {
        for (int j = n_cols - 1; j >= 0; j--) {
            double sum = 0;
            for (size_t k = 0; k < n_rows; k++) {
                sum += U[i * n_rows + k] * X[k * n_cols + j];
            }
            X[i * n_cols + j] = (y[i * n_cols + j] - sum) / U[i * n_rows + i];
        }
    }

    return X;
}

double* cholsky_solver(const double* L, const double* B, size_t n_rows, size_t n_cols){
    
    // Timer
    clock_t tik, tok;
    double time_taken;
    tik = clock();

    double* L_T= (double*) calloc(n_rows*n_cols,sizeof(double));

    transpose(n_rows, n_cols, L, L_T);
    
    double* Y = forward_substitusion(L,B,n_rows, n_cols);

    double* X = backwards_substitution(L_T,Y,n_rows, n_cols);
    free(Y);

    tok = clock();
    time_taken = (double)(tok - tik) / CLOCKS_PER_SEC;
    printf("C timer, Time taken for cholsky_solver: %f seconds\n", time_taken);
    
    return X;
}

对应的Python实现 #

仅使用Numpy的内置函数：

import numpy as np

# 假设A是对称正定矩阵，B是右侧矩阵
L, lower = np.linalg.cho_factor(A)
X = np.linalg.cho_solve((L, lower), B)

优化建议 #

1. 修复前向/后向替换的无效计算 #

这是最影响性能的核心问题：

• 前向替换中，L是下三角矩阵，第i行的k只需遍历0到i-1（k >= i时L[i*n_rows +k]为0，无需计算）。修改forward_substitusion的内层循环：

  for (size_t k = 0; k < i; k++) { // 替换原k < n_rows
      sum += L[i*n_rows + k] * Y[k*n_rows + j];
  }
  

• 后向替换中，U=L^T是上三角矩阵，第i行的k只需遍历i+1到n_rows-1（k <= i时U[i*n_rows +k]为0）。修改backwards_substitution的内层循环：

  for (size_t k = i+1; k < n_rows; k++) { // 替换原k < n_rows
      sum += U[i * n_rows + k] * X[k * n_cols + j];
  }
  

2. 避免显式转置操作 #

求解L^T*x=y时，无需分配内存存储L_T，直接通过索引访问L的转置元素：L[j*n_cols + i]等价于L_T[i*n_cols + j]。这样能节省内存分配和转置的时间开销，修改cholsky_solver：

double* cholsky_solver(const double* L, const double* B, size_t n_rows, size_t n_cols){
    clock_t tik, tok;
    double time_taken;
    tik = clock();

    double* Y = forward_substitusion(L,B,n_rows, n_cols);
    // 直接传入L，在backwards_substitution中处理转置逻辑
    double* X = backwards_substitution(L,Y,n_rows, n_cols);
    free(Y);

    tok = clock();
    time_taken = (double)(tok - tik) / CLOCKS_PER_SEC;
    printf("C timer, Time taken for cholsky_solver: %f seconds\n", time_taken);
    
    return X;
}

同时修改backwards_substitution，直接使用L的转置索引：

double* backwards_substitution(const double* L, const double* y, size_t n_rows, size_t n_cols) {
    double* X = (double*)calloc(n_rows * n_cols, sizeof(double));

    for (int i = n_rows - 1; i >= 0; i--) {
        for (int j = n_cols - 1; j >= 0; j--) {
            double sum = 0;
            for (size_t k = i+1; k < n_rows; k++) {
                // L[k*n_cols + i] 就是L^T[i*n_cols +k]
                sum += L[k * n_cols + i] * X[k * n_cols + j];
            }
            X[i * n_cols + j] = (y[i * n_cols + j] - sum) / L[i * n_cols + i];
        }
    }
    return X;
}

3. 内存布局与缓存优化 #

• 调整循环顺序：在Cholesky分解的三重循环中，将k循环放在最内层，同时保证内存访问连续。比如修改循环顺序为j -> i -> k，让L[i*n_cols +k]和L[j*n_cols +k]的访问更符合缓存局部性。
• 仅存储下三角部分：由于L是下三角矩阵，无需存储上三角的0元素，可将L压缩存储为一维数组（大小为n*(n+1)/2），减少内存占用并提升缓存命中率。

4. 开启编译器优化 #

编译时启用最高级别的优化选项，比如GCC使用-O3，Clang使用-O3 -march=native。这些选项会让编译器自动做循环展开、向量优化、寄存器分配等优化，大幅提升代码性能。

5. 调用BLAS/LAPACK库 #

Numpy的cho_factor和cho_solve底层调用了优化过的BLAS/LAPACK接口（如dpotrf做Cholesky分解，dpotrs做三角求解），这些库经过了高度优化，利用了CPU的SIMD指令、缓存优化甚至多线程。直接调用这些接口比手动实现高效得多，示例代码（使用OpenBLAS）：

#include <cblas.h>
#include <lapacke.h>

// Cholesky分解
lapack_int cholesky_factor_lapack(double* A, size_t n) {
    return LAPACKE_dpotrf(LAPACK_ROW_MAJOR, 'L', n, A, n);
}

// 求解
void cho_solve_lapack(double* L, double* B, size_t n, size_t nrhs) {
    LAPACKE_dpotrs(LAPACK_ROW_MAJOR, 'L', n, nrhs, L, n, B, n);
}

6. 减少内存分配开销 #

• 避免在函数内频繁使用calloc，可提前分配好内存并传入函数，或者使用栈内存（当矩阵不大时）。
• calloc会初始化内存为0，若不需要初始化（比如前向替换中Y初始值确实需要0）可保留，否则用malloc更快。

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