读 ncnn 源码（Ⅹ）：Winograd F(2×2,3×3) 的内核变换、选路与打包（含对比 packed sse）

TL;DR

• 入口：conv3x3s1_winograd23_transform_kernel() 把 3×3 卷积核先做 Winograd 核变换，再按 (M,K) 维度 分块打包到 AT，以便后续 batched-GEMM 高效消费。

• 选 tile：conv3x3s1_winograd_get_optimal_tile_mnk() 依据 L2 容量 + ISA 对齐 + 线程数，优先 不切 K，设定 (TILE_M,TILE_K[,TILE_N])。

• 打包：conv3x3s1_winograd_pack_A_tile() 将线程私有缓冲 A_tile 重新编织成 AT 的目标布局（以 batch=B=16 为内层 stride），方便微内核整块装载。

• 何时用 23/63：test_prefer_winograd{23,63} 给出基于 in/out/min(w,h) 的经验区间；3×3、s=1、d=1 且通道/空间合适时优先 Winograd。

• 对比：convolution_transform_kernel_packed_sse() 是非 Winograd路径的权重重排：按 elempack/out_elempack 将 (kw,kh,in,out) 变到向量友好布局。

  

---

1）Winograd 算法

这里给出Winograd F(2, 3) 算法的详细数学推导。这是理解 Winograd 算法最经典、最基础的案例。

目标：推导计算一维卷积 Y = conv(d, g) 的过程。

• 输出 Y：长度为 2，$Y = [y₀, y₁]$
• 卷积核 g：长度为 3，$g = [g₀, g₁, g₂]$
• 输入 d：长度为 $m+r-1 = 2+3-1 = 4$，$d = [d₀, d₁, d₂, d₃]$

直接卷积 (Direct Convolution)
根据卷积定义，我们有：

$y₀ = d₀g₀ + d₁g₁ + d₂g₂$

$y₁ = d₁g₀ + d₂g₁ + d₃g₂$

总共需要 6 次乘法 和 4 次加法。我们的目标就是减少乘法的次数。

---

第一步：将卷积问题转化为多项式乘法

这是 Winograd 算法的数学基石。我们可以将卷积核 g 和输入 d 分别表示为多项式：

• 卷积核多项式：$G(x) = g₀ + g₁x + g₂x²$
• 输入多项式：$D(x) = d₀ + d₁x + d₂x² + d₃x³$

将这两个多项式相乘，得到结果多项式 $Y(x) = G(x) \cdot D(x)$。展开这个乘法，我们会发现，结果多项式中 $x²$ 和 $x³$ 的系数 正好对应着卷积的输出 $y₀$ 和 $y₁$。（注意：根据不同的卷积定义和偏移，可能是不同次幂的系数，但其等价性是核心）。

为了简化推导，我们采用一种更通用的形式，定义输出多项式 $Y(x) = y₀ + y₁x$。我们的任务就是通过间接计算求得 $y₀$ 和 $y₁$。

---

第二步：选择求值点

为了求解一个次数为 $n$ 的多项式，我们需要在 $n+1$ 个点上对它进行求值。我们关心的多项式 $Y(x)$ 的最高次数为 $(m-1) + (r-1) = (2-1) + (3-1) = 3$。因此，我们需要在 $3+1=4$ 个点上进行求值。

Winograd 算法的精髓在于选择了这样一组求值点：$x = 0, 1, -1, \infty$。

• 选择 $0, 1, -1$ 是因为它们可以极大地简化计算。
• 选择 $\infty$ 是一个巧妙的技巧，它实际上代表着取多项式的最高次项系数。

---

第三步：推导变换矩阵

现在，我们在这四个点上，分别对输入多项式 $D(x)$ 和卷积核多项式 $G(x)$ 进行求值。

A. 推导输入变换矩阵 $B^T$

我们需要计算 $v = [D(0), D(1), D(-1), D(\infty)]^T$。

1. $D(0)$ = $d₀ + d₁(0) + d₂(0)² + d₃(0)³ = d₀$
2. $D(1)$ = $d₀ + d₁ + d₂ + d₃$
3. $D(-1)$ = $d₀ - d₁ + d₂ - d₃$
4. $D(\infty)$ = $d₃$ (取最高次项 $x³$ 的系数)

将这个过程写成矩阵形式：

$v = B^T \cdot d$

$\begin{bmatrix} D(0) \\ D(1) \\ D(-1) \\ D(\infty) \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 1 & 1 & 1 & 1 \\ 1 & -1 & 1 & -1 \\ 0 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} d₀ \\ d₁ \\ d₂ \\ d₃ \end{bmatrix}$

因此，我们得到了输入变换矩阵 $B^T$。
(注意：在 ncnn 等库的实际实现中，为了进一步优化加法，可能会使用一个稍有不同的 $B^T$ 矩阵，但其数学本质相同)

B. 推导卷积核变换矩阵 $G$

我们需要计算 $u = [G(0), G(1), G(-1), G(\infty)]^T$。

1. $G(0)$ = $g₀ + g₁(0) + g₂(0)² = g₀$
2. $G(1)$ = $g₀ + g₁ + g₂$
3. $G(-1)$ = $g₀ - g₁ + g₂$
4. $G(\infty)$ = $g₂$ (取最高次项 $x²$ 的系数)

写成矩阵形式：

$u = G \cdot g$

$\begin{bmatrix} G(0) \\ G(1) \\ G(-1) \\ G(\infty) \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ 1 & 1 & 1 \\ 1 & -1 & 1 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} g₀ \\ g₁ \\ g₂ \end{bmatrix}$

这就是卷积核变换矩阵 $G$。

C. 推导输出变换矩阵 $A^T$ (插值)

现在我们已经通过 $m = u \odot v$ 得到了结果多项式 $Y(x)$ 在四个点上的值：

• $m₀ = Y(0)$
• $m₁ = Y(1)$
• $m₂ = Y(-1)$
• $m₃ = Y(\infty)$

我们的目标是反解出 $y₀$ 和 $y₁$。标准的 F(2,3) 输出变换矩阵 $A^T$ 是通过拉格朗日插值法并结合特定的卷积定义反推出来的，其形式为：

$A^T = \begin{bmatrix} 1 & 1 & 1 & 0 \\ 0 & 1 & -1 & -1 \end{bmatrix}$

所以，输出的计算是：

$Y = A^T \cdot m$

$\begin{bmatrix} y₀ \\ y₁ \end{bmatrix} = \begin{bmatrix} 1 & 1 & 1 & 0 \\ 0 & 1 & -1 & -1 \end{bmatrix} \begin{bmatrix} m₀ \\ m₁ \\ m₂ \\ m₃ \end{bmatrix}$

即：

$y₀ = m₀ + m₁ + m₂$

$y₁ = m₁ - m₂ - m₃$

---

第四步：总结 Winograd F(2, 3) 完整流程

1. 离线步骤 (模型加载时)：

   • 用矩阵 $G$ 变换卷积核 $g$，得到 $u = G \cdot g$。

2. 在线步骤 (模型推理时)：

   • 输入变换：从输入 d 中提取 4 个元素，用矩阵 $B^T$ 变换，得到 $v = B^T \cdot d$。
   • 逐元素乘积：计算 $m = u \odot v$。这一步只包含 4 次乘法。
   • 输出变换：用矩阵 $A^T$ 变换 $m$，得到最终结果 $Y = A^T \cdot m$。

性能对比：

• 直接卷积：6 次乘法
• Winograd：4 次乘法
• 乘法运算量：减少了 (6-4)/6 ≈ 33.3%。

这个推导过程展示了 Winograd 算法如何通过将问题转换到多项式域，并利用精心选择的求值点，将变换的代价降至最低，从而显著减少核心计算中的乘法次数。对于二维情况 F(2x2, 3x3)，这个原理被应用两次，使得乘法次数从 36 次 减少到 16 次，带来了 2.25 倍的理论加速。

2）入口：分块 + 线程缓冲 + 目标布局

static void conv3x3s1_winograd23_transform_kernel(const Mat& kernel, Mat& AT, int inch, int outch, const Option& opt)
{
    const int M = outch;      // 输出通道
    const int K = inch;       // 输入通道
    const int B = 16;         // F(2,3) 产生 4×4 变换核 => batch = 16

    int TILE_M, TILE_N, TILE_K;
    conv3x3s1_winograd_get_optimal_tile_mnk(M, 0, K, B, TILE_M, TILE_N, TILE_K, opt.num_threads);

    Mat A_tileX(B * TILE_M * TILE_K, 1, opt.num_threads);  // 每线程私有缓冲
    AT.create(TILE_K * TILE_M, B, (K+TILE_K-1)/TILE_K, (M+TILE_M-1)/TILE_M);
    // 直观：AT 的 4D 视图 ~ [M块][K块][B(=16)][TILE_K*TILE_M] 的重排

并行外层按 M-block 切分，每个线程循环 K-block：

1. conv3x3s1_winograd23_transform_kernel_tile：把该 (M×K) 小块里的 每个 3×3 核变成 4×4；顺序写入线程缓冲 A_tile；
2. conv3x3s1_winograd_pack_A_tile：将 A_tile 打包到 AT 对应 tile 位置，使 batch（B=16）为最快变动维、内存连续，贴合 GEMM 微内核的装载习惯。

---

2）F(2×2,3×3) 的核变换：$g=G g G^{\mathsf T}$

// G = [[1,0,0],
//      [1/2, 1/2, 1/2],
//      [1/2,-1/2, 1/2],
//      [0,0,1]]

static inline void conv3x3s1_winograd23_transform_kernel_tile(
    const Mat& kernel, Mat& A, int inch, int i, int max_ii, int k, int max_kk)
{
    float* ptmp = A;              // 输出到线程私有 tile
    for (int ii = 0; ii < max_ii; ii++)
    {
        for (int kk = 0; kk < max_kk; kk++)
        {
            // 取出一个 3×3 核：outch = i+ii, inch = k+kk
            const float* k0 = (const float*)kernel + (i+ii)*inch*9 + (k+kk)*9;

            float tmp[4][3];      // 先做左乘 G：得到 4×3
            for (int m = 0; m < 3; m++)
            {
                float r0 = k0[0], r1 = k0[1], r2 = k0[2];
                tmp[0][m] = r0;
                tmp[1][m] = 0.5f*(r0 + r1 + r2);
                tmp[2][m = 0.5f*(r0 - r1 + r2);
                tmp[3][m] = r2;
                k0 += 3;          // 下一列
            }

            // 再右乘 G^T：对 tmp 的每行与 G^T 相乘 => 4×4
            for (int m = 0; m < 4; m++)
            {
                float r0 = tmp[m][0], r1 = tmp[m][1], r2 = tmp[m][2];
                float z0 = r0;
                float z1 = 0.5f*(r0 + r1 + r2);
                float z2 = 0.5f*(r0 - r1 + r2);
                float z3 = r2;

                ptmp[0]=z0; ptmp[1]=z1; ptmp[2]=z2; ptmp[3]=z3;
                ptmp += 4;        // 连续写 4 个，行优先
            }
        }
    }
}

要点

• 这段代码就是把 $G$ 和 $G^{\mathsf T}$ 的矩阵乘 手写展开：减少循环与访存，固定常数 0.5 便于编译器矢量化。
• 每个 3×3 核产出 4×4（共 16 个）“频域”权重；B=16 正是源于此。
• k0 += 3 的增量对应列访问（3×3 的行内 stride=3）。

---

3）把线程 tile 打包成 AT：batch 为内层

static void conv3x3s1_winograd_pack_A_tile(const Mat& A, Mat& AT, int batch, int max_ii, int max_kk)
{
    // A 中每个 3×3 核已变 4×4，共 batch=16 个系数 => Ktile * 16 是“列”跨度
    const int N = max_kk * batch;   // 每行的跨度

    for (int b = 0; b < batch; b++) // 以 b 为内层维度重排
    {
        float* pp = AT.row(b);
        int ii = 0;

        // NEON / AArch64 有更宽的块（8/4），这里省略……
        for (; ii + 1 < max_ii; ii += 2)         // 写两个 out channel
        {
            const float* p0 = (const float*)A + ii * N + b;

            for (int kk = 0; kk < max_kk; kk++)  // 穿过 Ktile
            {
                pp[0] = p0[0];
                pp[1] = p0[N];                   // 下一个 out channel（ii+1）
                p0 += batch;                     // 跳到下一组 (b 对应的) 变换系数
                pp += 2;
            }
        }
        // 余数分支：ii+0
        for (; ii < max_ii; ii++)
        {
            const float* p0 = (const float*)A + ii * N + b;
            for (int kk = 0; kk < max_kk; kk++)
            {
                *pp++ = *p0;
                p0 += batch;
            }
        }
    }
}

图像化理解

• 线程临时 A 的布局更“计算友好”（ii 外、kk×b 内）；
• AT 作为全局权重缓存，需要把 B=16 放到最内层（row(b)），这样连续装载就能拿到同一 (i,k) 的 16 个变换系数，喂给微内核做 batched-GEMM。
• NEON 的 8/4 宽块分支就是为了宽矢量一次写更多通道（ii 维度展开），减少循环控制与访存指令。

---

4）Winograd 的 tile 选择：不切 K，M 乘线程，兼顾 ISA 对齐

static void conv3x3s1_winograd_get_optimal_tile_mnk(int M,int N,int K,int B,int& TM,int& TN,int& TK,int nT)
{
    const int l2 = get_cpu_level2_cache_size()/sizeof(float);

    // 先解 K：尽量不切 K，根据 L2 容量和 ISA（aarch64/NEON/标量）
    // align 到 8/4/2 的倍数
    // …… 细节略

    // 再解 M：aarch64 基准 8、NEON 4、标量 2，然后 × 线程数，再根据 M 反向收敛对齐
    // …… 细节略

    // 若 N>0 再解 N（本核变换阶段 N=0，不涉及）
}

原则

• 优先不切 K：让每线程的工作集（TILE_M*TILE_K*B）尽量驻留 L2，减少跨核/跨层级访存。
• M 方向并行：TILE_M 基于 ISA 的向量宽度给一个下界（8/4/2），再放大到适配线程数，保证每个线程吃到一整块 M。
• 对齐：所有 tile 都会收敛到 8/4/2 的倍数，确保后续矢量转置/装载是整齐的。

---

5）何时偏好 Winograd（23 / 63）

test_prefer_winograd{23,63} 是 经验表（作者在特定平台 profile 得到），按 num_input/num_output/min(w,h) 给区间。粗略理解：

• F(2,3)（winograd23）适合 更大的空间/通道，提升更稳（代价是 B=16 的打包与内存）。
• F(6,3)（winograd63）在更小空间/通道时不一定划算，但在一些“中等”区域能大幅减少乘加次数。
• 一旦 stride != 1、dilation != 1 或 核不是 3×3，通常就不走 Winograd 了（回落到 sgemm/packed）。

---

6）对比：非 Winograd 的权重打包

static void convolution_transform_kernel_packed_sse(
    const Mat& w, Mat& w_tm, int in, int out, int kw, int kh, int ep, int oep)
{
    const int maxk = kw*kh;
    Mat w_r2 = w.reshape(maxk, in, out);
    w_tm.create(maxk, in/ep, out/oep, (size_t)4u*ep*oep, ep*oep);

    for (int q=0; q+oep-1<out; q+=oep)
    {
        float* g = w_tm.channel(q/oep);
        for (int p=0; p+ep-1<in; p+=ep)
        {
            for (int k=0; k<maxk; k++)
                for (int i=0; i<ep; i++)
                    for (int j=0; j<oep; j++)
                    {
                        const float* k00 = w_r2.channel(q+j).row(p+i);
                        *g++ = k00[k];
                    }
        }
    }
}

差异点

• 这里不做 $G g G^{\mathsf T}$ 的数学变换，只是把 (kw,kh,in,out) 的权重排列成 ep×oep 的向量块，方便 AVX/SSE 的打包内核消费。
• Winograd 路线的 AT 则是在此基础上多了一层 4×4 的“频域”维度（B=16），以及以 B 为内层的布局策略。

---

7）常见坑 & 实战建议

• F(2,3) 只适用 3×3 / s=1 / d=1：其他配置别强上。
• 通道与对齐：确保 outch、inch 有足够规模并能命中 ISA 对齐（8/4），否则 Winograd 的打包/转置开销可能吃掉收益。
• L2 容量敏感：自定义平台上最好 profile 一下 (TILE_M, TILE_K) 的选择。
• 内存形状：理解 AT.create(TK*TM, B, K/TK, M/TM) 的 4D 组织，有助于你写自定义微内核时一次性整块装载 B=16 的系数。

---

8) 对比

在 Winograd 算法的标准表示法 F(m, r) 中：

• r 是固定的卷积核尺寸，这里 r=3 (代表 3x3)。
• m 是输出tile的尺寸，也就是该算法一次计算能产生多大一块输出。

因此：

• 23 对应 F(2x2, 3x3)：一次计算产生 2x2 的输出块。
• 43 对应 F(4x4, 3x3)：一次计算产生 4x4 的输出块。
• 63 对应 F(6x6, 3x3)：一次计算产生 6x6 的输出块。

这个 m 的选择，直接决定了算法的性能特征。

特性 / 变体	F(2x2, 3x3) (23)	F(4x4, 3x3) (43)	F(6x6, 3x3) (63)
输出瓦块尺寸 (m)	2x2	4x4	6x6
输入瓦块尺寸 (m+r-1)	4x4	6x6	8x8
变换域尺寸	4x4	6x6	8x8
核心乘法次数	16	36	64
理论乘法加速比	36 / 16 = 2.25x	144 / 36 = 4.0x	324 / 64 = 5.06x
算术强度 (Arithmetic Intensity)	低	中	高
变换开销 (加法/访存)	最低	中等	最高
数值稳定性	最好	良好	可能稍差
最佳适用场景	小型特征图 (如网络末端)	通用，大多数情况下的甜点	大型特征图 (如网络前端)

---

9）一句话收尾

Winograd 在 ncnn 里的实现哲学是：先把 3×3 核变成 4×4 的“更便于乘加”的系数块，再把这些块打包成微内核喜欢的内存顺序——数学和工程两头都拿下，才有稳定的端上性能。

该封面图片由Ralf Ruppert在Pixabay上发布