在计算机图形学、物理模拟、机器人技术以及众多科学计算领域，3x3 和 4x4 矩阵的乘法运算极为常见且对性能至关重要。虽然现代 CPU 速度飞快，但在需要执行数百万次这类运算的场景下（例如实时渲染的每一帧），即使是微小的优化也能带来显著的性能提升。本报告将探讨几种加速这两种特定尺寸矩阵乘法的实用技术。

1. 基准：朴素算法 (Naive Algorithm)

我们首先回顾标准的矩阵乘法定义。对于两个矩阵 $A$ ($m \times n$) 和 $B$ ($n \times p$)，它们的乘积 $C = A \times B$ 是一个 $m \times p$ 的矩阵，其中每个元素 $C_{ij}$ 由下式给出：

$$ C_{ij} = \sum_{k=0}^{n-1} A_{ik} B_{kj} $$

对于 $3 \times 3$ 矩阵 ($m=n=p=3$) 或 $4 \times 4$ 矩阵 ($m=n=p=4$)，这对应于一个三重嵌套循环：

1
2
3
4
5
6
7
8
9
10
11
12
// 朴素 4x4 矩阵乘法示例 (C = A * B)
mat4 multiply_naive(const mat4& a, const mat4& b) {
    mat4 c = {}; // 初始化为零
    for (int i = 0; i < 4; ++i) {      // Row index for A and C
        for (int j = 0; j < 4; ++j) {  // Column index for B and C
            for (int k = 0; k < 4; ++k) { // Inner dimension index
                c.m[i][j] += a.m[i][k] * b.m[k][j];
            }
        }
    }
    return c;
}

该算法的时间复杂度为 $O(n^3)$。对于固定的 $n=3$ 或 $n=4$，这本身并不算糟糕，但其循环结构引入了不可忽视的开销（计数器增量、条件判断、分支预测），并且可能不利于现代 CPU 的缓存和流水线执行。

2. 优化技术一：循环展开 (Loop Unrolling)

循环展开是一种编译器优化技术，也可以手动实现，旨在通过减少或消除循环控制指令来降低开销。对于固定且较小的循环次数，我们可以完全展开循环。

2.1 完全展开 (3x3 矩阵)

对于 $3 \times 3$ 矩阵，乘法涉及计算 9 个结果元素。每个元素需要 3 次乘法和 2 次加法。总计 27 次乘法和 18 次加法。由于计算量固定且不大，我们可以完全展开所有循环，直接计算每个结果元素：

计算 $C_{00}$:
$$ C_{00} = A_{00}B_{00} + A_{01}B_{10} + A_{02}B_{20} $$
计算 $C_{01}$:
$$ C_{01} = A_{00}B_{01} + A_{01}B_{11} + A_{02}B_{21} $$
… 以此类推，直到 $C_{22}$。

代码原理:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
mat3 multiply_3x3_unrolled(const mat3& a, const mat3& b) {
    mat3 result; // C = A * B

    // 直接计算每个元素，无需循环
    // Row 0
    result.m[0][0] = a.m[0][0] * b.m[0][0] + a.m[0][1] * b.m[1][0] + a.m[0][2] * b.m[2][0];
    result.m[0][1] = a.m[0][0] * b.m[0][1] + a.m[0][1] * b.m[1][1] + a.m[0][2] * b.m[2][1];
    result.m[0][2] = a.m[0][0] * b.m[0][2] + a.m[0][1] * b.m[1][2] + a.m[0][2] * b.m[2][2];

    // Row 1 (类似计算)
    result.m[1][0] = a.m[1][0] * b.m[0][0] + a.m[1][1] * b.m[1][0] + a.m[1][2] * b.m[2][0];
    // ... result.m[1][1], result.m[1][2] ...

    // Row 2 (类似计算)
    result.m[2][0] = a.m[2][0] * b.m[0][0] + a.m[2][1] * b.m[1][0] + a.m[2][2] * b.m[2][0];
    // ... result.m[2][1], result.m[2][2] ...

    return result;
}

这种方法消除了所有循环开销，指令流水线可能更流畅。缺点是代码体积增大，但对于 3x3 来说通常是值得的，且具有良好的可移植性。

2.2 部分展开 (4x4 矩阵)

对于 $4 \times 4$ 矩阵，完全展开（计算 16 个元素，每个需要 4 次乘法和 3 次加法，总计 64 次乘法和 48 次加法）虽然可行，但代码会非常冗长。更常见的做法是仅展开最内层的循环（k 循环），它对应于点积计算。

代码原理:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
mat4 multiply_4x4_inner_unrolled(const mat4& a, const mat4& b) {
    mat4 result = {}; // Initialize to zero

    for (int i = 0; i < 4; ++i) {
        for (int j = 0; j < 4; ++j) {
            // 内层循环 k 被展开
            result.m[i][j] = a.m[i][0] * b.m[0][j] +
                             a.m[i][1] * b.m[1][j] +
                             a.m[i][2] * b.m[2][j] +
                             a.m[i][3] * b.m[3][j];
        }
    }
    return result;
}

这保留了外两层循环，但消除了最频繁执行的内层循环的开销，是一个很好的性能与代码复杂度之间的折衷。

3. 优化技术二：SIMD 向量化 (Vectorization)

现代 CPU 普遍支持 SIMD（Single Instruction, Multiple Data）指令集，如 SSE (Streaming SIMD Extensions) 和 AVX (Advanced Vector Extensions)。这些指令允许 CPU 在单个时钟周期内对多个数据元素（通常是 4 个 32 位浮点数或整数）执行相同的操作。

3.1 SIMD 与 4x4 矩阵

$4 \times 4$ 矩阵与 4 元素向量操作（如 SSE 的 128 位寄存器，可容纳 4 个 float）天然契合。我们可以将矩阵的行或列视为向量进行并行计算。

一种常见的 SIMD 策略是计算结果矩阵 $C$ 的一行 $C_i$：
$$ C_i = \sum_{k=0}^{3} A_{ik} B_k $$
其中 $A_{ik}$ 是标量（矩阵 $A$ 的元素），$B_k$ 是矩阵 $B$ 的第 $k$ 行（视为向量）。在 SIMD 实现中，这通常转化为：

  1. 将 $B$ 的 4 行加载到 4 个 SIMD 寄存器中。
  2. 对于 $A$ 的第 $i$ 行，将其元素 $A_{i0}, A_{i1}, A_{i2}, A_{i3}$ 逐个“广播”（复制）到 SIMD 寄存器的所有通道。
  3. 使用广播后的 $A_{ik}$ 值与加载的 $B_k$ 行向量进行并行乘法。
  4. 将 4 次乘法的结果累加起来，得到结果行 $C_i$。

代码原理 (SSE 示例):

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
#include <immintrin.h> // Or <xmmintrin.h>, <emmintrin.h> etc.

mat4 multiply_4x4_sse(const mat4& a, const mat4& b) {
    mat4 result;
    const float* pA = &a.m[0][0];
    const float* pB = &b.m[0][0];
    float* pResult = &result.m[0][0];

    for (int i = 0; i < 4; ++i) { // Calculate row i of result
        // Load rows of B into SSE registers (__m128 holds 4 floats)
        __m128 b_row0 = _mm_loadu_ps(&pB[0 * 4]); // Unaligned load row 0
        __m128 b_row1 = _mm_loadu_ps(&pB[1 * 4]); // Unaligned load row 1
        __m128 b_row2 = _mm_loadu_ps(&pB[2 * 4]); // Unaligned load row 2
        __m128 b_row3 = _mm_loadu_ps(&pB[3 * 4]); // Unaligned load row 3

        // Pointer to current row of A
        const float* a_row_ptr = &pA[i * 4];

        // Accumulator for the result row C[i]
        __m128 result_row;

        // C[i] = A[i][0] * B[0]
        result_row = _mm_mul_ps(_mm_set1_ps(a_row_ptr[0]), b_row0);
        // C[i] += A[i][1] * B[1]
        result_row = _mm_add_ps(result_row, _mm_mul_ps(_mm_set1_ps(a_row_ptr[1]), b_row1));
        // C[i] += A[i][2] * B[2]
        result_row = _mm_add_ps(result_row, _mm_mul_ps(_mm_set1_ps(a_row_ptr[2]), b_row2));
        // C[i] += A[i][3] * B[3]
        result_row = _mm_add_ps(result_row, _mm_mul_ps(_mm_set1_ps(a_row_ptr[3]), b_row3));

        // Store the calculated row into the result matrix
        _mm_storeu_ps(&pResult[i * 4], result_row); // Unaligned store
    }
    return result;
}

(_mm_set1_ps 用于广播标量，_mm_loadu_ps / _mm_storeu_ps 用于非对齐内存访问，_mm_mul_ps 和 _mm_add_ps 执行并行的乘法和加法)。

SIMD 实现通常能提供最高的性能，但代价是代码可移植性差（依赖特定指令集）和复杂性增加。

3.2 SIMD 与 3x3 矩阵

将 SIMD 应用于 $3 \times 3$ 矩阵比较棘手，因为 SIMD 寄存器通常是 4 宽的。需要进行数据填充、掩码操作或复杂的重排（shuffling），这可能引入额外的开销，抵消并行计算的优势。因此，对于 3x3 矩阵，完全循环展开通常是更实用、更高效的选择。

4. 编译器优化与标志

除了手动优化代码，编译器的优化能力也至关重要。

  • 优化级别: 务必启用优化标志，如 GCC/Clang 的 -O2 或 -O3，MSVC 的 /O2。这些标志会启用包括自动循环展开、指令重排、自动向量化（有时）在内的多种优化。
  • 目标架构: 使用 -march=native (GCC/Clang) 或 /arch:AVX2 (MSVC) 等标志，可以让编译器生成针对特定 CPU 指令集的优化代码，充分利用可用的 SIMD 功能，即使是在看似普通的 C++ 代码（如循环展开版本）上，编译器也可能生成 SIMD 指令。

5. 为何不用 Strassen 等算法？

Strassen 算法及其变种具有优于 $O(n^3)$ 的渐近复杂度（Strassen 约为 $O(n^{2.81})$）。然而，这些算法的常数因子和管理开销（递归、额外的加减法、内存分配）非常高。对于 $n=3$ 或 $n=4$ 这样的小尺寸，其开销远超其理论优势，实际性能通常不如经过优化的朴素算法。

6. 结论与建议

加速 3x3 和 4x4 矩阵乘法没有银弹，需要根据目标平台、性能需求和可接受的复杂性来选择策略：

  • 3x3 矩阵: 完全循环展开 通常是最佳选择，它提供了良好的性能提升，且代码相对直接，可移植性好。
  • 4x4 矩阵:
    如果追求极致性能且目标平台确定（例如 PC 游戏开发），手动 SIMD (SSE/AVX) 实现 是首选。

如果需要更好的可移植性或希望简化代码，内层循环展开 是一个可靠且有效的优化。
无论哪种手动优化，配合编译器的 -O2/-O3 和 -march=native//arch:AVX2 标志 都至关重要，以发挥硬件和编译器的全部潜力。

最后，任何性能优化都应基于实际测量 (Profiling)。在目标环境和典型负载下测试不同实现的性能，才能确定哪种方法真正“最快”。

License

本文采用 署名-非商业性使用-相同方式共享 4.0 国际 许可协议，转载请注明出处。