早前《Winograd卷积原理 | Hey~YaHei!)》一文已经介绍过Winograd的卷积原理，但有些细节似乎尚不明确——

• Winograd为何能够奏效？
• 它背后的数学原理是什么？
• 变换矩阵如何得到？

本文将借鉴《源于《孙子算经》的Cudnn | 知乎, Silicon Valley》和《wincnn/2464-supp.pdf | github, andravin》，深入剖析Winograd背后的数学原理。

中国同余定理

Winograd的基础来源于《孙子算经》中的中国同余定理（Chinese Remainder Theorem），同余定理是数论的重要内容（没学过数论的我，要读懂这个数学原理真的有些头疼QAQ）。

《孙子算经》

《孙子算经》有一章：

今有物不知其数，三三数之剩二，五五数之剩三，七七数之剩二，问物几何？
答曰：二十三。
术曰：三三数之剩二，置一百四十；五五数之剩三，置六十三，七七数之剩二，置三十，并之。得二百三十三，以二百一十减之，即得。凡三三数之剩一，则置七十；五五数之剩一，则置二十一；七七数之剩一，则置十五；一百六以上以一百五减之即得。

改用数学语言描述一下：

问题：已知同余方程组
$$\begin{cases} x \equiv 2(mod\ 3) \\ x \equiv 3(mod\ 5) \\ x \equiv 2(mod\ 7) \end{cases}$$
求x？

其中，$x \equiv a(mod\ m)$称为x与a模m同余，当$a < m$时，可以认为x模m余a；

解答：
$m_1 = 2 \times 70 = 140$
$m_2 = 3 \times 21 = 63$
$m_3 = 2 \times 15 = 30$
$x = (m_1 + m_2 + m_3)\ mod\ 105 = 23$

看到这，没读过数论的我感到十分震惊——为啥拿着余数乘几个数，最后加起来模一模怎么就可以推出被模的数呢？？

定理内容

这里不加证明地直接引出中国同余定理的内容：

假设整数$m_1, m_2, …, m_n$两两互素，则对于任意的整数$a_1, a_2, …, a_n$，同余方程组
$$\begin{cases} x \equiv a_1(mod\ m_1) \\ x \equiv a_2(mod\ m_2) \\ … \\ x \equiv a_n(mod\ m_n) \end{cases}$$
都存在整数解，且若$X, Y$都满足方程组，则必有$x \equiv Y (mod\ N)$，其中$N=\prod^n_{i=1}m_i$
$$x \equiv \sum^n_{i=1} a_i \times \frac{N}{m_i} \times [(\frac{N}{m_i})^{-1}]_{m_i} (mod\ N)$$
其中$[(\frac{N}{m_i})^{-1}]_{m_i}$是一个整数，使得$\frac{N}{m_i} \times [(\frac{N}{m_i})^{-1}]_{m_i} \equiv 1(mod\ m_i)$

对《孙子算经》中的例子来说，就是
$$
m_1 = 3, m_2 = 5, m_3 = 7 \\
a_1 = 2, a_2 = 3, a_3 = 2 \\
N = m_1 \times m_2 \times m_3 = 105 \\
\frac{N}{m_1} = 35, \frac{N}{m_2} = 21, \frac{N}{m_3} = 15 \\
[(\frac{N}{m_1})^{-1}]_{m_1} = 2, [(\frac{N}{m_2})^{-1}]_{m_2} = 1, [(\frac{N}{m_3})^{-1}]_{m_3} = 1
$$

原始的定理只针对整数，事实上，该定理同样适用于多项式。
接下来我们关注一下定理里的两个细节：

如何证明互素？

互素，又称互质。若两个整数的公约数只有1，那么称它们为互素（质）整数。
当然这个定义可以引申到多项式上，比如$x$、$x+1$、$x-1$、$x^2+1$就是典型的互素多项式。

证明互素的核心是证明它们的公约数只有1，问题可以转化成最大公约数为1，说起最大公约数的计算，那就不得不说欧几里得算法啦！

欧几里得算法俗称辗转相除法，用大数除以小数得到余数（这里商是无意义的），然后迭代计算，用上一轮的除数、余数分别作为新一轮的被除数和除数，直到最后能够整除，最后一轮的除数即为原始的两个数的最大公约数。
以1997和615为例：
$1997\ mod\ 615 = 152$
$615\ mod\ 152 = 7$
$152\ mod\ 7 = 5$
$7\ mod\ 5 = 2$
$5\ mod\ 2 = 1$
$2\ mod\ 1 = 0$
故1997与615的最大公约数为1，记为gcd(1997, 615) = 1

至于，多项式之间如何做除法，可以直接用小学的长除法，除到最后余数的最高次幂小于除数的最高次幂即可。

以$x^4 + x^3 - x^2 + 3x + 2$除以$x^3 + 3x^2 + 3x + 2$为例：
（latex没法写长除公式，这里我只能用分数的形式来表述）
$$\begin{align} \frac{x^4 + x^3 - x^2 + 3x + 2}{x^3 + 3x^2 + 3x + 2} &= \frac{-2x^3 - 4x^2 + x + 2}{x^3 + 3x^2 + 3x + 2} + (x) \\ &= \frac{2x^2 + 7x + 6}{x^3 + 3x^2 + 3x + 2} + (x-2) \end{align}$$
那么这里$x^4 + x^3 - x^2 + 3x + 2$除以$x^3 + 3x^2 + 3x + 2$，商为$(x-2)$，余数为$(2x^2 + 7x + 6)$

如何求解系数$[(\frac{N}{m_i})^{-1}]_{m_i}$？

根据定理，系数$[(\frac{N}{m_i})^{-1}]_{m_i}$满足$\frac{N}{m_i} \times [(\frac{N}{m_i})^{-1}]_{m_i} \equiv 1(mod\ m_i)$，我们假设
$$\left( \frac{N}{m_i} \times [(\frac{N}{m_i})^{-1}]_{m_i} \right) \div m_i = y \cdots 1$$
那么有
$$\frac{N}{m_i} \times [(\frac{N}{m_i})^{-1}]_{m_i} - m_i y = 1$$
同时注意到，$\frac{N}{m_i}$与$m_i$必定互素，则有
$$\frac{N}{m_i} \times [(\frac{N}{m_i})^{-1}]_{m_i} - m_i y = 1 = gcd(\frac{N}{m_i}, m_i)$$
构成了一个形如$ax + by = gcd(a, b)$的贝祖公式，按照贝祖定理，可用扩展欧几里得算法求解出$x$和$y$，也即$[(\frac{N}{m_i})^{-1}]_{m_i}$和$y$（当然y没啥意义，我们要的是那个系数！）。

扩展欧几里得算法除了能和欧几里得算法一样求得最大公约数之外，还能得到模逆元；
其实非常简单，与递归实现的欧几里得算法相比，无非是每次递归都逐步恢复出逆元x跟y
以python实现为例：

# def euclid(a, b):
#     while b != 0:
#         a, b = b, a % b
#     return a

def euclid(a, b):
    if b == 0:
        return a
    else:
        return euclid(b, a % b)

def ext_euclid(a, b):
    if b == 0:
        return 1, 0, a
    else:
        x, y, q = ext_euclid(b, a % b)
        x, y = y, (x - (a // b) * y)
        return x, y, q

Winograd算法

前边中国同余定理已经解释明白了，接下来我们试着理解一下Winograd的数学原理（以$F(2,3)$为例）。
假设两个离散序列，
$$
\begin{align}
d[n] &= [d_0, d_1] \\
g[n] &= [g_0, g_1, g_2]
\end{align}
$$
经过Z变换变成多项式形式（为了方便起见，这里写作$x$而非$Z$），
$$
\begin{align}
d(x) &= d_1x + d_0 \\
g(x) &= g_2x^2 + g_1x + g_0
\end{align}
$$
此时，$d[n] * g[n] = Z^{-1}(d(x) \cdot g(x))$（时域卷积 => Z域的乘法），记
$$
\begin{align}
y &= d(x) \cdot g(x) \\
&= (d_1x + d_0)(g_2x^2 + g_1x + g_0)
\end{align}
$$
于是两个序列的卷积变成了多项式的乘法（当然，乘积也是个多项式），我们继续往多项式版本的中国同余定理上套。
选取K个互素的一次多项式$m_i(x)$，使得$m(x) = \prod^{K-1}_{i=0}{m_i(x)}$的阶数等于$d(x) \cdot g(x)$，那就可以构造出同余方程
$$y(x) - R_{m(x)}[y(x)] \equiv g(x)d(x)\ (mod\ m(x))$$
其中，$R_{m(x)}[y(x)]$是$y(x)$除以$m(x)$的余数。
在这个例子中，至少需要3个互素一次多项式，不妨取
$$
\begin{align}
m_0(x) &= x \\
m_1(x) &= x + 1 \\
m_2(x) &= x - 1
\end{align}
$$
那么有$m(x) = \prod{m_i(x)} = x(x+1)(x-1) = x^3 - x$，
除此之外，我们再记一个特殊的$m_3(x) = x - \infty$，没有具体的意义，该多项式是为了凑足m的阶数，使得$m’(x) = m_3(x) m(x) = x(x+1)(x-1)(x-\infty)$阶数高于$d(x) \cdot g(x)$，从而可以构造出同余方程
$$y(x) \equiv g(x)d(x)\ (mod\ m’(x))$$
记$M_i(x) = \frac{m(x)}{m_i(x)}$，则有
$$
\begin{align}
M_0(x) &= -1 + x^2 \\
M_1(x) &= -x + x^2 \\
M_2(x) &= x + x^2 \\
M_3(x) &= -x + x^3
\end{align}
$$
依次抽取$x^0, x^1, x^2, x^3$的系数组成行向量并堆叠构造出变换矩阵
$$
B =
\begin{bmatrix}
-1 & 0 & 0 & 0 \\
0 & -1 & 1 & -1 \\
1 & 1 & 1 & 0 \\
0 & 0 & 0 & 1
\end{bmatrix}
$$
注意
$$y(x)\ mod\ M_3(x) = y(x)\ mod\ m(x)$$

接下来考虑一个重要的同余方程性质：

若$a \equiv b(mod\ m_i), i \in \{1, 2, …, n\}$
则$a \equiv b(mod\ [m_1, m_2, …, m_n])$，
其中$[m_1, m_2, …, m_n]$表示$m_1, m_2, …, m_n$的最小公倍数

利用该性质（因为$m_1, m_2, …, m_n$两两互素，所以$m’(x)$就是它们的最小公倍数），又根据同余方程的乘法规则

$$a \equiv b(mod\ m), c \equiv d(mod\ m) => ac \equiv bd(mod\ m)$$

原同余方程$y(x) \equiv g(x)d(x)\ (mod\ m’(x))$可以转化为同余方程组
$$
\begin{cases}
y(x) \equiv g(x)d(x)\ (mod\ m_0(x)) \\
y(x) \equiv g(x)d(x)\ (mod\ m_1(x)) \\
y(x) \equiv g(x)d(x)\ (mod\ m_2(x)) \\
y(x) \equiv g(x)d(x)\ (mod\ m_3(x))
\end{cases}
$$
的求解。原方程组又可以展开成
$$
\begin{cases}
y(x) \equiv [g(x)\ mod\ m_0(x)] \cdot [d(x)\ mod\ m_0(x)]\ (mod\ m_0(x)) \\
y(x) \equiv [g(x)\ mod\ m_1(x)] \cdot [d(x)\ mod\ m_1(x)]\ (mod\ m_1(x)) \\
y(x) \equiv [g(x)\ mod\ m_2(x)] \cdot [d(x)\ mod\ m_2(x)]\ (mod\ m_2(x)) \\
y(x) \equiv [g(x)\ mod\ m_3(x)] \cdot [d(x)\ mod\ m_3(x)]\ (mod\ m_3(x)) \\
\end{cases}
$$
噢？好像变成了g(x)和d(x)，橘势似乎明朗了起来=。=
$$
\begin{align}
&\begin{cases}
g^{(0)}(x) = g(x)\ mod\ m_0(x) = g_0 \\
g^{(1)}(x) = g(x)\ mod\ m_1(x) = g_0 - g_1 + g_2 \\
g^{(2)}(x) = g(x)\ mod\ m_2(x) = g_0 + g_1 + g_2 \\
g^{(3)}(x) = g(x)\ mod\ m_3(x) = g_2
\end{cases} \\
&\begin{cases}
d^{(0)}(x) = d(x)\ mod\ m_0(x) = d_0 \\
d^{(1)}(x) = d(x)\ mod\ m_1(x) = d_0 - d_1 \\
d^{(2)}(x) = d(x)\ mod\ m_2(x) = d_0 + d_1 \\
d^{(3)}(x) = d(x)\ mod\ m_3(x) = d_1 \\
\end{cases}
\end{align}
$$
依次抽取$d_0, d_1$的系数组成列向量并堆叠构造出变换矩阵
$$
A =
\begin{bmatrix}
1 & 0 \\
1 & -1 \\
1 & 1 \\
0 & 1
\end{bmatrix}
$$
那么回过头看看中国同余定理，
$$x \equiv \sum^n_{i=1} a_i \times \frac{N}{m_i} \times [(\frac{N}{m_i})^{-1}]_{m_i} (mod\ N)$$
在这个例子里表示为
$$y(x) \equiv \sum^3_{i=0} g^{(i)}(x)d^{(i)}(x) \times M_i(x) \times N_i(x)\ (mod\ m’(x))$$
显然我们还差一个$N_i(x)$，根据贝祖公式$m_i(x)n_i(x) + M_i(x)N_i(x) = 1$用扩展欧几里得算法依次求得：
注意：此处《wincnn/2464-supp.pdf | github, andravin》给出的贝祖公式有误！！
$$
\begin{cases}
N_0(x) = -1 \\
N_1(x) = \frac{1}{2} \\
N_2(x) = \frac{1}{2} \\
N_3(x) = 1
\end{cases}
$$
而
$$
\begin{cases}
N_0(x)g^{(0)}(x) = -g_0 \\
N_1(x)g^{(1)}(x) = \frac{1}{2}g_0 - \frac{1}{2}g_1 + \frac{1}{2}g_2 \\
N_2(x)g^{(2)}(x) = \frac{1}{2}g_0 + \frac{1}{2}g_1 + \frac{1}{2}g_2 \\
N_3(x)g^{(3)}(x) = g_2
\end{cases}
$$
抽取$g_0, g_1, g_2$的系数组成列向量并堆叠构造出矩阵
$$
G =
\begin{bmatrix}
-1 & 0 & 0 \\
\frac{1}{2} & -\frac{1}{2} & \frac{1}{2} \\
\frac{1}{2} & \frac{1}{2} & \frac{1}{2} \\
0 & 0 & 1
\end{bmatrix}
$$

我们代回去验证一下，
$$
\begin{cases}
g^{(0)}(x)d^{(0)}(x) \times M_0(x) \times N_0(x) = -g_0d_0x^2 + g_0d_0 \\
g^{(1)}(x)d^{(1)}(x) \times M_1(x) \times N_1(x) = \frac{1}{2}(g_0-g_1+g_2)(d_0-d_1)x^2 - \frac{1}{2}(g_0-g_1+g_2)(d_0-d_1)x \\
g^{(2)}(x)d^{(2)}(x) \times M_2(x) \times N_2(x) = \frac{1}{2}(g_0+g_1+g_2)(d_0+d_1)x^2 + \frac{1}{2}(g_0+g_1+g_2)(d_0+d_1)x \\
g^{(3)}(x)d^{(3)}(x) \times M_3(x) \times N_3(x) = g_2d_1x^3 - g_2d_1x
\end{cases}
$$
代入求得
$$y(x) \equiv g_2d_1x^3 + (g_2d_0 + g_1d_1)x^2 + (g_0d_1 + g_1d_0)x + g_0d_0\ (mod\ m’(x))$$
$m’(x)$阶数为4，大于3，则有
$$y(x) = g_2d_1x^3 + (g_2d_0 + g_1d_1)x^2 + (g_0d_1 + g_1d_0)x + g_0d_0 = g(x)d(x)$$

同样的由于$m’(x)$阶数较高，原同余方程可以推出
$$y(x) = \sum^3_{i=0} g^{(i)}(x)d^{(i)}(x) \times M_i(x) \times N_i(x)$$
改写成矩阵形式并做Z反变换变成
$$
\begin{bmatrix} y_0 \\ y_1 \\ y_2 \\ y_3 \end{bmatrix} =
\left(
\begin{bmatrix}
-1 & 0 & 0 \\
\frac{1}{2} & -\frac{1}{2} & \frac{1}{2} \\
\frac{1}{2} & \frac{1}{2} & \frac{1}{2} \\
0 & 0 & 1
\end{bmatrix}
\begin{bmatrix} g_0 \\ g_1 \\ g_2 \end{bmatrix}
\right)
\odot
\left(
\begin{bmatrix}
-1 & 0 & 0 & 0 \\
0 & -1 & 1 & -1 \\
1 & 1 & 1 & 0 \\
0 & 0 & 0 & 1
\end{bmatrix}
\begin{bmatrix}
1 & 0 \\
1 & -1 \\
1 & 1 \\
0 & 1
\end{bmatrix}
\begin{bmatrix} d_0 \\ d_1 \end{bmatrix}
\right)
$$
也即
$$y= (Gg) \odot (BAd)$$
……

啊！卡住了，我已经推导不下去……

1. 上述Winograd的推导针对的是离散序列的卷积，我该如何推广到深度学习的卷积（即互相关函数）上去？
2. 为何m输入r参数的离散序列卷积推导出来的变换矩阵G、A、B可以应用在m输出r参数的互相关函数上，最后变成$y= A^T[(Gg) \odot (B^Td)]$的？