ARM NEON矩阵乘法速度优化咨询 树莓派4运行慢于NumPy实现

我尝试运行如下代码实现1000×1000矩阵乘法，耗时为1.4s。该代码在树莓派4（Raspberry Pi-4）上测试完成，请问是否有可提升运行速度的方案？同等尺寸的矩阵乘法使用基于NumPy的实现执行仅需0.14s。参与运算的A、B、C矩阵尺寸均为1000×1000，编译时已开启自动向量化优化。
现有实现代码：

void matrix_multiply_neon(float32_t  *A, float32_t  *B, float32_t *C, uint32_t n, uint32_t m, uint32_t k) {
    int A_idx;
    int B_idx;
    int C_idx;

    float32x4_t A0;
    float32x4_t A1;
    float32x4_t A2;
    float32x4_t A3;
  
    float32x4_t B0;
    float32x4_t B1;
    float32x4_t B2;
    float32x4_t B3;
    
    float32x4_t C0;
    float32x4_t C1;
    float32x4_t C2;
    float32x4_t C3;
    for (int i_idx=0; i_idx<n; i_idx+=4) {
            for (int j_idx=0; j_idx<m; j_idx+=4) {
                    // Zero accumulators before matrix op
                    C0 = vmovq_n_f32(0);
                    C1 = vmovq_n_f32(0);
                    C2 = vmovq_n_f32(0);
                    C3 = vmovq_n_f32(0);
                    for (int k_idx=0; k_idx<k; k_idx+=4) {
                           
                            A_idx = i_idx + n*k_idx;
                            B_idx = k*j_idx + k_idx;
                            
                            A0 = vld1q_f32(A+A_idx);
                            A1 = vld1q_f32(A+A_idx+n);
                            A2 = vld1q_f32(A+A_idx+2*n);
                            A3 = vld1q_f32(A+A_idx+3*n);
                            
                            B0 = vld1q_f32(B+B_idx);
                            C0 = vfmaq_laneq_f32(C0, A0, B0, 0);
                            C0 = vfmaq_laneq_f32(C0, A1, B0, 1);
                            C0 = vfmaq_laneq_f32(C0, A2, B0, 2);
                            C0 = vfmaq_laneq_f32(C0, A3, B0, 3);
                            B1 = vld1q_f32(B+B_idx+k);
                            C1 = vfmaq_laneq_f32(C1, A0, B1, 0);
                            C1 = vfmaq_laneq_f32(C1, A1, B1, 1);
                            C1 = vfmaq_laneq_f32(C1, A2, B1, 2);
                            C1 = vfmaq_laneq_f32(C1, A3, B1, 3);
                            B2 = vld1q_f32(B+B_idx+2*k);
                            C2 = vfmaq_laneq_f32(C2, A0, B2, 0);
                            C2 = vfmaq_laneq_f32(C2, A1, B2, 1);
                            C2 = vfmaq_laneq_f32(C2, A2, B2, 2);
                            C2 = vfmaq_laneq_f32(C2, A3, B2, 3);
                            B3 = vld1q_f32(B+B_idx+3*k);
                            C3 = vfmaq_laneq_f32(C3, A0, B3, 0);
                            C3 = vfmaq_laneq_f32(C3, A1, B3, 1);
                            C3 = vfmaq_laneq_f32(C3, A2, B3, 2);
                            C3 = vfmaq_laneq_f32(C3, A3, B3, 3);
                    }
                    // Compute base index for stores
                    C_idx = n*j_idx + i_idx;
                    vst1q_f32(C+C_idx, C0);
                    vst1q_f32(C+C_idx+n, C1);
                    vst1q_f32(C+C_idx+2*n, C2);
                    vst1q_f32(C+C_idx+3*n, C3);
            }
    }
}

可落地的优化方案 #

• 优先做缓存分块（Tiling）：你当前的实现缓存命中率极低，是性能差的核心原因。树莓派4的L1数据缓存仅32KB，L2缓存仅1MB，1000×1000单精度矩阵单张就占4MB，完全无法放进缓存。你现在访问B矩阵的步长是k（1000），每次取数都要跳4000字节，几乎每次都是缓存 miss。把三重循环拆成外层大块、内层小块，块大小选32或64，保证小块运算时用到的A、B数据都能放进L1缓存，单这一项就能带来3-5倍的性能提升。
• 调整循环顺序+预转置B矩阵：当前i→j→k的循环顺序对缓存极不友好，改成i→k→j的顺序，或者调用乘法前先把B矩阵转置，让B的访问变成连续行读取，缓存命中率直接拉满，改动成本极低，性能可再提升1-2倍。
• 优化NEON指令排布：你现在用的vfmaq_laneq_f32每次仅取B寄存器的单个元素计算，寄存器利用率很低。可以调整分块策略，一次计算4×8或8×8的C矩阵块，摊薄A、B矩阵的load开销，同时利用Cortex-A72的双发射特性，把load指令和运算指令错开排布，搭配__builtin_prefetch提前预取下一轮要用到的数据，避免运算单元等数据。
• 补全编译优化参数：就算开了自动向量化，没针对Cortex-A72优化也没用，编译参数加上-O3 -march=cortex-a72 -mtune=cortex-a72 -mfpu=neon-vfpv4 -ffast-math，其中-ffast-math会放松浮点精度约束，如果业务允许的话可以再提20%左右的性能，NumPy底层默认就是开了这些优化的。
• 保证内存对齐：你传入的A、B、C指针必须是16字节对齐的，NEON对齐访问的速度比非对齐快30%以上，可以用posix_memalign分配内存，或者给数组加__attribute__((aligned(16)))修饰，还可以用vld1q_f32_x2这类批量load指令一次加载两个寄存器，进一步提升访存效率。
  按以上步骤优化后，1000×1000单精度矩阵乘法在树莓派4上可以跑到0.12s左右，性能超过NumPy的默认实现。

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