稀疏矩阵乘法算法优化咨询（优先C++实现）

核心优化思路 #

针对你遇到的矩阵B列访问开销大的问题，核心优化方向是调整矩阵B的存储格式，让列访问变成高效的行访问，同时结合缓存局部性优化：

1. 转置矩阵B为CSR格式 #

将原本行存储（CSR）的B转置为B^{T（同样用CSR存储），这样原问题中需要获取的B(k,j)等价于B}T(j,k)。此时当遍历A的行i的非零元素A(i,k)时，可以直接访问B^T的行k的非零元素，通过行索引快速定位，避免遍历B的行查找列j的元素。

2. 改用CSC格式存储矩阵B #

CSC（列压缩存储）是专门针对列访问优化的格式：

• 用col_ptr数组记录每个列的非零元素在values和row_indices数组中的起始位置
• 用row_indices存储非零元素的行号，values存储对应数值
  这样要获取B的列j的所有非零元素，直接通过col_ptr[j]和col_ptr[j+1]就能定位到对应的元素范围，无需遍历。

3. 分块优化（利用缓存局部性） #

将大矩阵拆分为固定大小的子块（比如32x32，匹配CPU缓存线大小），优先计算子块内的乘法，减少缓存失效次数，提升整体运算效率。

C++代码示例 #

1. 定义CSR和CSC数据结构 #

#include <vector>
#include <unordered_map>
#include <algorithm>

// CSR格式：行压缩存储
struct CSRMatrix {
    int rows, cols;
    std::vector<double> values;
    std::vector<int> row_ptr; // row_ptr[i]是第i行第一个非零元素在values中的索引，row_ptr[rows]是总非零元素数
    std::vector<int> col_indices;
};

// CSC格式：列压缩存储
struct CSCMatrix {
    int rows, cols;
    std::vector<double> values;
    std::vector<int> col_ptr; // col_ptr[j]是第j列第一个非零元素在values中的索引，col_ptr[cols]是总非零元素数
    std::vector<int> row_indices;
};

2. 将CSR矩阵转置为CSR格式（得到B^T） #

CSRMatrix transposeCSR(const CSRMatrix& mat) {
    CSRMatrix transposed;
    transposed.rows = mat.cols;
    transposed.cols = mat.rows;
    transposed.values.resize(mat.values.size());
    transposed.col_indices.resize(mat.col_indices.size());
    transposed.row_ptr.resize(mat.cols + 1, 0);

    // 统计每列的非零元素数量
    for (int idx = 0; idx < mat.values.size(); ++idx) {
        int col = mat.col_indices[idx];
        transposed.row_ptr[col + 1]++;
    }

    // 计算前缀和得到row_ptr
    for (int j = 0; j < mat.cols; ++j) {
        transposed.row_ptr[j + 1] += transposed.row_ptr[j];
    }

    // 填充转置矩阵的元素
    std::vector<int> temp_ptr = transposed.row_ptr;
    for (int i = 0; i < mat.rows; ++i) {
        for (int idx = mat.row_ptr[i]; idx < mat.row_ptr[i + 1]; ++idx) {
            int col = mat.col_indices[idx];
            int pos = temp_ptr[col]++;
            transposed.values[pos] = mat.values[idx];
            transposed.col_indices[pos] = i;
        }
    }

    return transposed;
}

3. 基于B转置的优化乘法实现 #

CSRMatrix multiplyCSRWithTransposedB(const CSRMatrix& A, const CSRMatrix& B) {
    CSRMatrix B_T = transposeCSR(B);
    CSRMatrix C;
    C.rows = A.rows;
    C.cols = B.cols;
    C.row_ptr.resize(A.rows + 1, 0);

    // 先统计每行的非零元素数量
    std::vector<std::unordered_map<int, double>> row_temp(A.rows);
    for (int i = 0; i < A.rows; ++i) {
        for (int idxA = A.row_ptr[i]; idxA < A.row_ptr[i + 1]; ++idxA) {
            int k = A.col_indices[idxA];
            double valA = A.values[idxA];
            // 访问B^T的行k（对应原B的列k）
            for (int idxBT = B_T.row_ptr[k]; idxBT < B_T.row_ptr[k + 1]; ++idxBT) {
                int j = B_T.col_indices[idxBT];
                double valB = B_T.values[idxBT];
                row_temp[i][j] += valA * valB;
            }
        }
        C.row_ptr[i + 1] = C.row_ptr[i] + row_temp[i].size();
    }

    // 填充结果矩阵的values和col_indices
    C.values.reserve(C.row_ptr.back());
    C.col_indices.reserve(C.row_ptr.back());
    for (int i = 0; i < A.rows; ++i) {
        for (const auto& pair : row_temp[i]) {
            C.col_indices.push_back(pair.first);
            C.values.push_back(pair.second);
        }
    }

    return C;
}

4. 基于CSC格式B的乘法实现 #

// 将CSR转换为CSC
CSCMatrix CSRtoCSC(const CSRMatrix& mat) {
    CSCMatrix csc;
    csc.rows = mat.rows;
    csc.cols = mat.cols;
    csc.values.resize(mat.values.size());
    csc.row_indices.resize(mat.col_indices.size());
    csc.col_ptr.resize(mat.cols + 1, 0);

    for (int idx = 0; idx < mat.values.size(); ++idx) {
        int col = mat.col_indices[idx];
        csc.col_ptr[col + 1]++;
    }

    for (int j = 0; j < mat.cols; ++j) {
        csc.col_ptr[j + 1] += csc.col_ptr[j];
    }

    std::vector<int> temp_ptr = csc.col_ptr;
    for (int i = 0; i < mat.rows; ++i) {
        for (int idx = mat.row_ptr[i]; idx < mat.row_ptr[i + 1]; ++idx) {
            int col = mat.col_indices[idx];
            int pos = temp_ptr[col]++;
            csc.values[pos] = mat.values[idx];
            csc.row_indices[pos] = i;
        }
    }

    return csc;
}

CSRMatrix multiplyCSRwithCSCB(const CSRMatrix& A, const CSCMatrix& B) {
    CSRMatrix C;
    C.rows = A.rows;
    C.cols = B.cols;
    C.row_ptr.resize(A.rows + 1, 0);

    std::vector<std::unordered_map<int, double>> row_temp(A.rows);
    for (int i = 0; i < A.rows; ++i) {
        for (int idxA = A.row_ptr[i]; idxA < A.row_ptr[i + 1]; ++idxA) {
            int k = A.col_indices[idxA];
            double valA = A.values[idxA];
            // 直接访问B的列k的所有非零元素
            for (int idxB = B.col_ptr[k]; idxB < B.col_ptr[k + 1]; ++idxB) {
                int j = B.row_indices[idxB];
                double valB = B.values[idxB];
                row_temp[i][j] += valA * valB;
            }
        }
        C.row_ptr[i + 1] = C.row_ptr[i] + row_temp[i].size();
    }

    C.values.reserve(C.row_ptr.back());
    C.col_indices.reserve(C.row_ptr.back());
    for (int i = 0; i < A.rows; ++i) {
        for (const auto& pair : row_temp[i]) {
            C.col_indices.push_back(pair.first);
            C.values.push_back(pair.second);
        }
    }

    return C;
}

性能对比说明 #

• 转置B的方案：预处理开销是O(nnz(B))，但后续乘法中每个元素访问都是O(1)定位，适合多次复用B的场景。
• CSC存储B的方案：无需转置，直接针对列访问优化，单次乘法的效率更高，适合单次乘法场景。
• 分块优化可以叠加到上述方案中，进一步提升缓存命中率，对于超大稀疏矩阵效果显著。

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