深度学习（十二）——Winograd（2）-CSDN博客

本文链接：https://blog.csdn.net/antkillerfarm/article/details/78829889

最大公约数和Euclidean algorithm（续）

Euclidean algorithm的步骤如下图所示：

这里写图片描述

1.假设 $a>b$ ，则令 $c:=a \mod{b}$ 。

2.如果 $c=0$ ，则 $GCD(a,b)=b$ 。

3.否则令 $a:=b,b:=c$ ，并返回到第1步。

这个算法应该是Euclid记述的前人成果，因为更早的Eudoxus of Cnidus曾提到过这个算法。

Eudoxus of Cnidus，公元前390年～公元前337年，古希腊几何学家、天文学家和地理学家。柏拉图同时代最杰出的数学家。《几何原本》卷Ⅴ和卷Ⅻ主要来自欧多克索斯的工作。

然而，小学课本不使用Euclidean algorithm是有原因的，除了Euclidean algorithm本身相对复杂之外，短除法能同时搞定最大公约数和最小公倍数（Least common multiple），这也是它的教学优势所在。

Euclidean algorithm作为最古老的算法之一，被收录进Knuth的巨著TAOCP。这里的算法，指的是那些根据一定的规则来一步步执行的运算。

Bézout’s identity和Extended Euclidean algorithm

Euclidean division还可以表示为如下形式：

如果 $a'=a \mod(b)$ ，则 $a=\lfloor\frac{a}{b}\rfloor b+a'$ 。

这里的 $\lfloor x \rfloor$ 是向下取整的意思。

我们可以把上述定义扩展到负数。例如：

- 103 = - 2 \times 60 + 17 \Rightarrow (- 103) mod 60 = 17

$-103=-2\times 60+17 \Rightarrow (-103) \mod{60}=17$

注意：余数永远 $\ge 0$ 。

还可以把GCD的定义扩展为： $GCD(a,0)=a$ ，即任何整数都能整除0。

Euclidean division的定义扩展之后，则有Bézout’s identity。

Bézout’s identity：若a,b是非0整数， $d=GCD(a,b)$ ，则存在整数x,y使得 $ax+by=d$ 。证明略。

例如：

m 0 = 3, M 0 = 20, (- 1) 20 + 7 (3) = 1

$m_0 = 3 , M_0 = 20 , ( − 1 ) 20 + 7 ( 3 ) = 1$

Étienne Bézout，1730~1783，法国数学家。法国科学院院士。

至于如何求解x,y，这就要用到Extended Euclidean algorithm了。

首先，我们考虑边界情况，当 $b=0$ 时，原方程可化为 $ax=GCD(a,0)$ ，则：

{x = 1 y = 0 (1)

$\begin{cases} x=1 \\ y=0 \\ \end{cases}\tag{1}$

接着，设 $a'=b,b'=a \mod{b}$ ，则由Euclidean algorithm可得： $GCD(a,b)=GCD(b,a \mod b)$ ，因此：

a ​' ​ ​ x ​' ​ ​ + b ​' ​ ​ y ​' ​ ​ = G C D (a ​' ​ ​, b ​' ​ ​) = G C D (b, a mod b)

$a′x′+b′y′=GCD(a′,b′)=GCD(b,a \mod{b})$

a x + b y = a ​' ​ ​ x ​' ​ ​ + b ​' ​ ​ y ​' ​ ​ = b x ​' ​ ​ + (a mod b) y ​' ​ ​ = b x ​' ​ ​ + (a - ⌊ a b ⌋ b) y ​' ​ ​

$ax+by=a′x′+b′y′=bx′+(a\mod{b})y′=bx′+(a−\lfloor\frac{a}{b}\rfloor b)y′$

整理得到：

a x + b y = a y ​' ​ ​ + b (x ​' ​ ​ - ⌊ a b ⌋ y ​' ​ ​)

$ax+by=ay′+b(x′−\lfloor\frac{a}{b}\rfloor y′)$

对比系数，可得：

{x = y' y = x' - ⌊ a b ⌋ y ​' ​ ​ (2)

$\begin{cases} x=y' \\ y=x'−\lfloor\frac{a}{b}\rfloor y′ \\ \end{cases}\tag{2}$

公式1和2合到一起，就是一种迭代算法，也就是Extended Euclidean algorithm了。

从上面的讨论可知，Extended Euclidean algorithm实际上只在Euclidean algorithm之上前进了很小的一步，它的主要内核还是来源于Euclid。

但Euclid之所以不能更进一步，则主要是受制于负数的概念。虽然现在的小学高年级课本中，已经引入了负数，古代中国、印度、阿拉伯也很早就用到了负数，但是西方差不多要到文艺复兴时期，才逐渐接受了负数的概念。

不过反例也是有的，比如无理数，其它文明貌似根本就没有关注过它和有理数究竟有何区别…

参考：

https://blog.sengxian.com/algorithms/gcd-extgcd

欧几里德算法与扩展欧几里德算法

中国剩余定理

Chinese remainder theorem算是初等数论中，一个非常重要的定理了。（初等数论意指使用不超过高中程度的初等代数处理的数论问题，其最主要的工具包括整数的整除性与同余。）

CRT最早出自中国四世纪成书的古书《孙子算经》。著名的娱乐圈学霸关晓彤同学所攻克的“鸡兔同笼问题”，就出自该书。

CRT的内容为：

设 $m_i$ 为两两互质（pairwise coprime）的大于1的整数， $a_i$ 为任意整数，则存在x满足：

x \equiv a 1 ⋮ x \equiv a k (mod m 1) (mod m k) (1) (2) (3)

$\begin{align} x \equiv a_1 & \pmod{m_1} \\ \quad \vdots \\ x \equiv a_k &\pmod{m_k} \end{align}$

如果 $0\le x < M,M=\prod_{i=1}^k m_i$ ，则该x是唯一的。

CRT的存在性证明略。

这里以如下简单的例子，来讲讲如何求解x。

x x x \equiv 0 (mod 3) \equiv 3 (mod 4) \equiv 4 (mod 5) . (4) (5) (6)

$\begin{align} x &\equiv 0 \pmod{3} \\ x &\equiv 3 \pmod{4} \\ x &\equiv 4 \pmod{5}. \end{align}$

这个问题的穷举法需要遍历0到M的所有整数，这显然是个十分低效的算法。因此无论手算还是计算机算，基本都不用穷举法。

再来介绍一下筛法（Sieving）：

1.首先对 $m_i$ 按降序排序。

2.选择最大的模（这里为5）和对应的 $a_i$ （这里为4）。

{% highlight text %}
4 mod 4 → 0. Continue
4 + 5 = 9 mod 4 →1. Continue
9 + 5 = 14 mod 4 → 2. Continue
14 + 5 = 19 mod 4 → 3. OK, continue by considering remainders modulo 3 and adding 5×4 = 20 each time
19 mod 3 → 1. Continue
19 + 20 = 39 mod 3 → 0. OK, this is the result.
{% endhighlight %}

筛法对于M较小的情况，是非常高效的，因此手算一般都采用该法。但是，筛法的复杂度是指数级的，对于M较大的情况，并不好用。

CRT虽然只是初等数论的基本定理，但应用范围很广，Lagrange interpolation（一阶多项式插值）、Hermite interpolation（多阶多项式插值）和Dedekind’s theorem，都用到了CRT。

多项式的Euclidean division和GCD

我们可以仿照整数Euclidean division定义多项式的Euclidean division，如下面的竖式所示：

x 2 + 1 x + 3 x - 3) x 3 - 2 x 2 + 0 x - 4 ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ x 3 - 3 x 2 + 0 x - 4 - - - - - - - - - - - - - - - - - + x 2 + 0 x - 4 + x 2 - 3 x - 4 - - - - - - - - - - - - - + 3 x - 4 + 3 x - 9 - - - - - - - - + 5

$\begin{array}{r} x^2 + {\color{White}1}x + 3\\ x-3\overline{) x^3 - 2x^2 + 0x - 4}\\ \underline{x^3 - 3x^2 \color{White}{ + 0x - 4}}\\ +x^2 + 0x \color{White}{ - 4}\\ \underline{+x^2 - 3x \color{White}{ - 4}}\\ +3x - 4\\ \underline{+3x - 9}\\ +5 \end{array}$

上式也可写为横式：

x 3 - 2 x 2 - 4 = (x - 3) (x 2 + x + 3)              q (x) + 5    r (x)

${x^3 - 2x^2 - 4} = (x-3)\,\underbrace{(x^2 + x + 3)}_{q(x)} +\underbrace{5}_{r(x)}$

其中的 $r(x)$ 即为余数。

同样的可以定义多项式的GCD：

x 2 + 7 x + 6 = (x + 1) (x + 6)

$x^2 + 7x + 6 = (x + 1)(x + 6)$

x 2 - 5 x - 6 = (x + 1) (x - 6)

$x^2 − 5x − 6 = (x + 1)(x − 6)$

则两多项式的GCD为 $(x + 1)$ 。

多项式的CRT

CRT亦可改为如下等效形式：

c = (\sum i = 0 k c i N i M i) mod M

$c=\left(\sum_{i=0}^kc_iN_iM_i\right)\mod{M}$

其中 $m_i$ 两两互质， $c_i=R_{m_i}[c],M=\prod_{i=0}^km_i,M_i=M/m_i$ ， $N_i,n_i$ 是方程 $N_i M_i + n_i m_i = GCD ( M_i , m_i ) = 1$ 的解。

显然这里的 $N_i,n_i$ 可以使用Extended Euclidean algorithm求解。

例如：

m 0 = 3, M 0 = 20, (- 1) 20 + 7 (3) = 1

$m_0 = 3 , M_0 = 20 , ( − 1 ) 20 + 7 ( 3 ) = 1$

m 1 = 4, M 1 = 15, (- 1) 15 + (4) 4 = 1

$m_1 = 4 , M_1 = 15 , ( − 1 ) 15 + ( 4 ) 4 = 1$

m 2 = 5, M 2 = 12, (- 2) 12 + (5) 5 = 1

$m_2 = 5 , M_2 = 12 , ( − 2 ) 12 + ( 5 ) 5 = 1$

稍加扩展，可得到多项式版本的CRT：

c (p) = (\sum i = 0 k c (i) (p) N (i) (p) M (i) (p)) mod M (p) (3)

$c(p)=\left(\sum_{i=0}^kc^{(i)}(p)N^{(i)}(p)M^{(i)}(p)\right)\mod{M(p)}\tag{3}$

其中 $m^{(i)}(p)$ 两两互质， $c^{(i)}(p)=R_{m^{(i)}}[c(p)],M(p)=\prod_{i=0}^km^{(i)}(p),M^{(i)}(p)=M(p)/m^{(i)}(p)$ ， $N^{(i)}(p)$ 是方程 $N^{(i)}(p) M^{(i)}(p) + n^{(i)}(p) m^{(i)}(p) = GCD ( M^{(i)}(p) , m^{(i)}(p)) = 1$ 的解。

Winograd algorithm

下面以一个2x3的卷积为例，介绍一下Winograd algorithm的做法。

2x3卷积的多项式形式为：

h (p) = h 0 + h 1 p, x (p) = x 0 + x 1 p + x 2 p 2, s (p) = h (p) x (p)

$h(p)=h_0+h_1p,x(p)=x_0+x_1p+x_2p^2,s(p)=h(p)x(p)$

这里引入多项式（polynomial）的度（degree）的概念：多项式中包含的最高次项的次数，被称为多项式的度。

例如，上面的 $h(p)$ 的degree为1，而 $x(p)$ 的degree为2，而 $s(p)$ 的degree为3。

和Cook-Toom algorithm一样，Winograd algorithm也是一个构造式的算法。

Step 1：首先要构造一个degree大于等于3的多项式：

m (p) = m (0) (p) m (0) (p) \dots m (k) (p)

$m(p)=m^{(0)}(p)m^{(0)}(p)\cdots m^{(k)}(p)$

其中的 $m^{(i)}(p)$ 两两互质。

这里为了简单起见，不妨令 $m(p)=p(p-1)(p+1)$ ，并使用Extended Euclidean algorithm构建如下计算表格：

i	$m (i) (p)$ $m^{(i)}(p)$	$M (i) (p)$ $M^{(i)}(p)$	$n (i) (p)$ $n^{(i)}(p)$	$N (i) (p)$ $N^{(i)}(p)$
0	$p$ $p$	$p^{2} - 1$ $p^2-1$	$p$ $p$	$- 1$ $-1$
1	$p - 1$ $p-1$	$p 2 + p$ $p^2+p$	$- 1 2 (p + 2)$ $-\frac{1}{2}(p+2)$	$1 2$ $\frac{1}{2}$
2	$p + 1$ $p+1$	$p 2 - p$ $p^2-p$	$- 1 2 (p - 2)$ $-\frac{1}{2}(p-2)$	$1 2$ $\frac{1}{2}$

Step 2：使用如下公式计算 $h^{(i)}(p),x^{(i)}(p)$ ：

h (i) (p) = h (p) mod m (i) (p)

$h^{(i)}(p)=h(p)\mod{m^{(i)}(p)}$

x (i) (p) = x (p) mod m (i) (p)

$x^{(i)}(p)=x(p)\mod{m^{(i)}(p)}$

计算过程如下：

h (0) (p) = h 0, x (0) (p) = x 0

$h^{(0)}(p)=h_0,x^{(0)}(p)=x_0$

h (1) (p) = h 0 + h 1, x (1) (p) = x 0 + x 1 + x 2

$h^{(1)}(p)=h_0+h_1,x^{(1)}(p)=x_0+x_1+x_2$

h (2) (p) = h 0 - h 1, x (2) (p) = x 0 - x 1 + x 2

$h^{(2)}(p)=h_0-h_1,x^{(2)}(p)=x_0-x_1+x_2$

Step 3：使用如下公式计算 $s'^{(i)}(p)$ ：

s' (i) (p) = h (i) (p) x (i) (p) mod m (i) (p)

$s'^{(i)}(p)=h^{(i)}(p)x^{(i)}(p)\mod{m^{(i)}(p)}$

计算过程如下：

s' (0) (p) = h 0 x 0

$s'^{(0)}(p)=h_0x_0$

s' (1) (p) = (h 0 + h 1) (x 0 + x 1 + x 2)

$s'^{(1)}(p)=(h_0+h_1)(x_0+x_1+x_2)$

s' (2) (p) = (h 0 - h 1) (x 0 - x 1 + x 2)

$s'^{(2)}(p)=(h_0-h_1)(x_0-x_1+x_2)$

Step 4：根据公式3计算余数 $s'(p)$ ，并利用如下公式计算被除数 $s(p)$ ：

s (p) = s' (p) + h N - 1 x L - 1 m (p)

$s(p)=s'(p)+h_{N-1}x_{L-1}m(p)$

计算过程如下：

s (p) = s' (p) + h 1 x 2 m (p) = s' (0) (- p 2 + 1) + s ' ( 1 ) 2 (p 2 + p) + s ' ( 2 ) 2 (p 2 - p) + h 1 x 2 m (p 3 - p) = s' (0) + p (s ' ( 1 ) 2 - s ' ( 2 ) 2 - h 1 x 2) + p 2 (- s' (0) + s ' ( 1 ) 2 + s ' ( 2 ) 2) + p 3 (h 1 x 2) (13) (14) (15)

$\begin{align} s(p)&=s'(p)+h_1x_2m(p) \\ &=s'(0)(-p^2+1)+\frac{s'(1)}{2}(p^2+p)+\frac{s'(2)}{2}(p^2-p)+h_1x_2m(p^3-p)\\ &=s'(0)+p(\frac{s'(1)}{2}-\frac{s'(2)}{2}-h_1x_2)+p^2(-s'(0)+\frac{s'(1)}{2}+\frac{s'(2)}{2})+p^3(h_1x_2) \end{align}$

这里用4个乘法和7个加法，替代了6个乘法和2个加法。

总的来说，Winograd algorithm是一个很复杂的算法，但是结论却很简单。因此，在具体的IC实现中，一般只针对特定常用尺寸的kernel，应用相应的结论即可。