早前《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)$$
……
啊！卡住了，我已经推导不下去……
数学渣表示我可太难了>_<，推了四五天只能推出变换矩阵的推导过程和离散序列上卷积的Winograd算法。。至今不知道该怎么推广到深度学习的卷积（即互相关函数）的$y=A^T[(Gg)\odot(B^Td)]$上去
原文链接：《Winograd数学原理 | Hey~YaHei!》
