基于C语言优化矩阵计算：求矩阵与其转置相乘的高效方法及公式

Hey there! Let's dive into optimizing matrix-transpose multiplication in C—this is a problem where combining smart math with cache-aware code can make a huge difference. First, let's clarify the core operation we're targeting: given an m×n matrix A, we want to compute either B = A × Aᵀ (resulting in an m×m symmetric matrix) or C = Aᵀ × A (resulting in an n×n symmetric matrix). Let's break down the optimal mathematical ideas first, then tie them to C-specific optimizations.

1. Leverage the Symmetry of the Result #

The biggest mathematical win here is recognizing that both AAᵀ and AᵀA are symmetric matrices. That means for all indices i, j:

• B[i][j] = B[j][i] (for B = AAᵀ)
• C[i][j] = C[j][i] (for C = AᵀA)

Instead of calculating every element from scratch, you only need to compute the upper (or lower) triangle of the result matrix, then copy those values to their mirror positions. This cuts the number of arithmetic operations by roughly 50%—a massive efficiency boost for large matrices.

Mathematically, for B = AAᵀ, this looks like:

For i from 0 to m-1:
    For j from i to m-1:
        B[i][j] = sum_{k=0 to n-1} A[i][k] * A[j][k]
    For j from 0 to i-1:
        B[i][j] = B[j][i]

2. Block Matrix Decomposition (Cache Locality) #

Modern CPUs rely heavily on cache performance—accessing data from L1/L2 cache is orders of magnitude faster than main memory. Naive matrix multiplication leads to frequent cache misses when accessing non-contiguous elements.

The fix is to split your original matrix A into smaller, cache-sized sub-blocks. Mathematically, if we partition A into blocks of size b×b (where b is chosen to fit in your CPU's cache, e.g., 32 or 64), then:

• AAᵀ can be computed as the sum of products of these blocks and their transposes.

For example, if A is split into blocks A_pq (p ranges over row blocks, q over column blocks), then:

B_pr = sum_{q=0 to num_col_blocks-1} A_pq × A_rqᵀ

This way, each block fits entirely into cache, minimizing the number of times data has to be fetched from main memory. In C, this translates to reordering nested loops to iterate over blocks first, then elements within blocks.

3. Vectorized Inner Product Calculation #

Each element B[i][j] is the dot product of the i-th row and j-th row of A (for AAᵀ). Dot products are perfect candidates for SIMD (Single Instruction, Multiple Data) instructions, which let you compute multiple arithmetic operations in parallel.

Mathematically, the dot product can be broken into chunks that fit into SIMD registers. For example, with AVX-2 (256-bit registers), you can compute 8 single-precision floating-point multiplications at once, then accumulate the results.

In C, you can use compiler intrinsics (like _mm256_mul_ps and _mm256_add_ps for floats) to directly leverage these instructions. The core idea is to vectorize the inner sum:

sum = 0.0
for k in 0 to n-1 step 8:
    vec_a = load 8 elements from A[i][k..k+7]
    vec_b = load 8 elements from A[j][k..k+7]
    sum += dot_product(vec_a, vec_b)

• Loop Order: Always prioritize accessing data in row-major order (since C stores matrices row-wise). For block-based code, the order should be: block rows → block columns → inner block rows → inner block columns.
• Avoid Aliasing: Use restrict pointers to help the compiler optimize memory accesses (it tells the compiler that pointers don't overlap).
• Compiler Optimizations: Enable flags like -O3, -mavx2, or -mfma (depending on your CPU) to let the compiler auto-vectorize and optimize your code.
• Unroll Inner Loops: Manually unroll small inner loops (or let the compiler do it with -funroll-loops) to reduce loop overhead.

Here's a simplified C code skeleton that combines symmetry and blocking to boost performance:

#include <string.h>

void multiply_transpose(float *A, float *B, int m, int n, int block_size) {
    // Initialize result matrix to zero
    memset(B, 0, m * m * sizeof(float));

    for (int i_block = 0; i_block < m; i_block += block_size) {
        for (int j_block = i_block; j_block < m; j_block += block_size) {
            // Compute upper triangle block and mirror to lower triangle
            for (int i = i_block; i < i_block + block_size && i < m; i++) {
                for (int j = j_block; j < j_block + block_size && j < m; j++) {
                    float sum = 0.0f;
                    for (int k = 0; k < n; k++) {
                        sum += A[i * n + k] * A[j * n + k];
                    }
                    B[i * m + j] = sum;
                    if (i != j) {
                        B[j * m + i] = sum;
                    }
                }
            }
        }
    }
}

This skeleton cuts redundant calculations via symmetry and uses blocking to keep data in cache. For even better performance, you'd add SIMD intrinsics to the inner k loop.

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