Winograd Convolution Algorithm

• Winograd Convolution Algorithm is a family of fast algorithms for small fixed-size convolutions that reduce arithmetic operations using input, filter, and inverse transforms.
• It achieves significant computational savings, e.g., a 2.25× reduction for F(2,3) and up to 4× for F(4,3), by replacing direct multiply-accumulate operations with structured transforms.
• The method balances performance gains with challenges like numerical instability at larger tile sizes, mitigated by optimized point selection, mixed-precision, and novel techniques such as NOVA.

The Winograd Convolution Algorithm is a family of fast algorithms for implementing small fixed-size convolutions with provably minimal arithmetic complexity. Used extensively in both classical signal processing and deep learning, the method replaces direct multiply-accumulate (MAC) evaluation of convolutions with a sequence of input and filter transforms, elementwise products (Hadamard multiplications), and inverse transforms. This approach exploits the algebraic structure of finite impulse response (FIR) convolution and can achieve substantial performance gains on modern CPUs, GPUs, FPGAs, and AI inference ASICs.

1. Mathematical Foundations and Minimal Filtering

Winograd's minimal filtering algorithm for 1D convolution, denoted $F(m, r)$, computes $m$ outputs of the convolution between an input segment of length $m+r-1$ and a filter of length $r$ using only $m+r-1$ multiplications—achieving the theoretical lower bound for linear convolution. In matrix form, the computation is:

$y = A^T \bigl[ (Gg) \circ (B^T d) \bigr]$

where $g \in \mathbb{R}^r$ is the filter, $d \in \mathbb{R}^{m+r-1}$ is the data segment, $G\in\mathbb{R}^{(m+r-1)\times r}$, $B^T\in\mathbb{R}^{(m+r-1)\times(m+r-1)}$, $A^T\in\mathbb{R}^{m\times(m+r-1)}$, and $"\circ"$ denotes the Hadamard (elementwise) product.

In two dimensions, this methodology generalizes to produce $m \times m$ output tiles,

$Y = A^T\left[ (GgG^T) \circ (B^T d B) \right] A$

where, for a $r \times r$ kernel, the transform and inverse-transform matrices are typically derived via polynomial interpolation at $\omega = m + r - 1$ distinct points, forming (partial) Vandermonde matrices (Lavin et al., 2015, Tong et al., 2021, Liu et al., 2021).

Key Formulaic Dimensions:

• Direct 2D convolution on an $m \times m$ tile with an $r \times r$ kernel: $m^2 r^2$ multiplications.
• Winograd F(m, r): $(m + r - 1)^2$ multiplications per tile, corresponding to a reduction ratio of $m^2 r^2 / (m + r - 1)^2$; e.g., for $F(2, 3)$, $36/16 = 2.25 \times$ fewer multiplications (Liu et al., 2021, Lavin et al., 2015, Tong et al., 2021).

2. Construction of Transform Matrices

The transform matrices $B, G, A$ are derived by sampling polynomials at chosen interpolation points (the “point selection” problem). For a standard $F(m, r)$ configuration, $\omega = m + r - 1$ points are selected (e.g., $\{0, \pm1, \infty\}$ for $F(4,3)$), and the following are constructed:

• $B^T$: data transform (Vandermonde on input points)
• $G$: filter transform (scaled Vandermonde/Lagrange coefficients)
• $A^T$: output or inverse transform (inverse/interpolation matrix)

The explicit construction is via the method of Chinese Remainder Theorem for polynomials, or, equivalently, Lagrange/Toom–Cook interpolation (Lavin et al., 2015, Liu et al., 2021). Extensions to arbitrary $m, r$ require careful numerical choices of interpolation points; for example, (Barabasz et al., 2018, Alam et al., 2022, Lohia, 20 Dec 2025) detail how real-valued, rational, or even complex points can yield favorable conditioning and reduced error amplification in $G, B, A$.

3. Arithmetic Complexity, Speedup, and Limits

The principal advantage of the Winograd algorithm lies in operation-count reduction. For $3 \times 3$ convolution kernels, the commonly used $F(2,3)$ and $F(4,3)$ configurations achieve approximately $2.25\times$ and $4.0\times$ reduction in multiplies, respectively (Liu et al., 2021, Tong et al., 2021, Lavin et al., 2015). However, as the tile size increases (larger $m$), the benefit saturates due to transform overhead and exponentially growing numerical errors from ill-conditioned transform matrices (Vandermonde matrices exhibit poor conditioning for large $m$).

In practice, optimal performance is attained for small to medium $m$ (typically $m=2$ or $4$ in floating-point, $m=2$ in INT8), with larger tile sizes being compromised by floating-point and quantization inaccuracies (Barabasz et al., 2018, Lavin et al., 2015).

4. Numerical Stability, Quantization, and Conditioning

The convergence to minimal arithmetic complexity via Winograd comes at a cost: the conditioning of the transform matrices. Numerical instability is endemic for larger tile sizes or when using low-precision arithmetic (e.g., FP16, INT8). The worst-case error grows with the product of the norms and condition numbers of the $G, B, A$ transforms—often exponentially with tile size (Barabasz et al., 2018, Lohia, 20 Dec 2025, Alam et al., 2022). Several remedies have been developed:

• Optimized point selection: Using real-valued or symmetric points in structured forms (e.g., $\{\pm1/c, \pm c\}$) can drastically reduce error by simplifying the entries and the conditioning of $G, B, A$ (Alam et al., 2022).
• NOVA algorithm: Treats the point selection as a continuous optimization problem to discover fractional-valued points with dramatically improved conditioning; achieves up to $415\times$ condition number reduction for $F(8,3)$ in 1D (over $172,484\times$ in 2D) and stabilizes FP16/INT8 inference at large tile sizes (Lohia, 20 Dec 2025).
• Mixed-precision and pairwise summation: Using higher-precision for transforms, balanced summation trees (e.g., Huffman/priority pairing), and careful accumulation strategies can further reduce rounding error (Barabasz et al., 2018).
• Winograd-aware quantization and training: Embedding the quantization pipeline inside the Winograd transform/inverse, and optionally relaxing $G,B,A$ as trainable parameters (“flex” transforms) during training, has enabled robust INT8 inference—restoring accuracy losses up to 75 percentage points for otherwise unstable tiles (Fernandez-Marques et al., 2020, Barabasz, 2020).

5. Extensions: Large Kernels, Arbitrary Strides, and Nested Constructions

While the canonical algorithm is limited to small (e.g., $3\times3$) kernels and unit stride, several advances have generalized Winograd for broader applicability:

• Nested Winograd: For large kernel convolutions ($R\times R$, e.g., $7\times7$ up to $31\times31$), the nested Winograd decomposition recursively applies small-tab Winograd along each axis, reducing total multiplications by up to $10.5\times$ over linear decomposition—critical for accelerator throughput (Jiang et al., 2021).
• Decomposable Winograd Method (DWM): Decomposes large kernel and/or stride convolution into multiple small-shape Winograd-amenable sub-convolutions; achieves consistent $\sim2\times$ speedup with preserved accuracy and negligible error growth (Huang et al., 2020).
• Beyond Toom–Cook: Extensions allow use of higher-degree (nonlinear) irreducible polynomials (e.g., $a^2+1$), reducing error growth at the same multiply count, crucial for fp16/bf16 (Barabasz et al., 2019).
• Residue Number System (RNS): Enables low-precision integer domain Winograd by mapping all arithmetic into parallel RNS lanes, with CRT-based reconstruction—achieving up to $7.03\times$ theoretical reduction with zero loss in quantized CNNs (Liu et al., 2020).

6. Implementation: Hardware, Software, and Optimization

Efficient realization of Winograd convolution requires careful architecture-specific tiling, data reuse, and transform fusion:

• GPU and CPU: High-performance implementations fuse the entire Winograd pipeline (transforms, GEMM, inverse) into a single kernel/block, exploiting cache locality and vector units; memory layouts (e.g., “z-shape” packing, NHWC storage) are optimized for contiguous access (Gui et al., 2024, Lavin et al., 2015).
• FPGA/ASIC: Systolic arrays with parameterized Winograd PEs (WinoPEs), kernel sharing via unified selector bits, and BRAM-based streaming achieve DSP utilization up to $1.33$ GOPS/DSP and overall throughput over $3.1$ TOPS. Custom tap-wise quantization and power-of-two scaling, as well as fault-tolerance enhancements, are reported in state-of-the-art DSA designs (Liu et al., 2021, Andri et al., 2022, Xue et al., 2022).
• Sparse Winograd: Pruning in the Winograd domain leads to $>90\%$ sparsity in kernels, making possible $5.4\times$ speedup over dense direct convolution, especially impactful on general-purpose CPUs (Li et al., 2017).

7. Applications, Performance, and Future Directions

The Winograd algorithm is foundational in high-throughput, energy-efficient DNN accelerators, being the default kernel for $3\times3$ convolutions in cuDNN, OpenVINO, MKL-DNN, and numerous FPGA/ASIC toolflows. Layerwise and end-to-end throughput improvements of $1.5$–$4\times$ are empirically achieved, with energy reduction and resilience under both hardware faults and algorithmic quantization (Tong et al., 2021, Liu et al., 2021, Xue et al., 2022, Andri et al., 2022).

Ongoing research extends Winograd to large-kernel CNNs (nested decompositions), arbitrary architectures (nonlinear polynomials), energy/fault-aware systems, and extreme quantization (FP16/INT8/RNS), with automated search for optimal transform points (e.g., NOVA) (Lohia, 20 Dec 2025, Jiang et al., 2021, Alam et al., 2022). The method is increasingly vital for unlocking efficient, robust inference on next-generation AI hardware.

---

References: (Liu et al., 2021, Lavin et al., 2015, Tong et al., 2021, Jiang et al., 2021, Lohia, 20 Dec 2025, Barabasz et al., 2018, Alam et al., 2022, Huang et al., 2020, Andri et al., 2022, Li et al., 2017, Liu et al., 2020, Fernandez-Marques et al., 2020, Xue et al., 2022, Barabasz et al., 2019, Gui et al., 2024, Barabasz, 2020).