自動求導的二三事

----------------------------------------------------------------------------------------------------

每次搬到知乎也是夠麻煩的，不能直接貼markdown的代碼，更多博客見DarkScope從這裡開始

----------------------------------------------------------------------------------------------------

知乎上看到一個回答，說是自己學習神經網路的時候都是自己對公式求導，現在常見的DL庫都可以自動求導了。這個想必實現過神經網路的同學都有體會，因為神經網路的back-propagation演算法本質上就是求導鏈式法則的累積，所以學習這部分的時候就是推來推去，推導對了，那演算法你也就掌握了。

粗粗一想，只要能把所有操作用有向圖構建出來，通過遞歸去實現自動求導似乎很簡單，一時興起寫了一些代碼，整理成博客記錄一下。

[tips]完整代碼見這裡just.dark的代碼

#動手

首先我們需要一個基礎類，所有有向圖的節點都會有下面兩個方法`partialGradient()`是對傳入的變數求偏導，返回的同樣是一個圖。

`expression()`是用於將整個式子列印出來

class GBaseClass(): def __init__(self, name, value, type): self.name = name self.type = type self.value = value pass def partialGradient(self, partial): pass def expression(self): pass

從圖1可以看出來，我們主要有三種Class，常量Constant，變數Variable以及運算元Operation，最簡單的常量：

G_STATIC_VARIABLE = {} class GConstant(GBaseClass): def __init__(self, value): global G_STATIC_VARIABLE try: G_STATIC_VARIABLE[counter] += 1 except: G_STATIC_VARIABLE.setdefault(counter, 0); self.value = value self.name = "CONSTANT_" + str(G_STATIC_VARIABLE[counter]) self.type = "CONSTANT" def partialGradient(self, partial): return GConstant(0) def expression(self): return str(self.value)

可以看到，我們為常量設置了自增的name，只需要傳入value即可定義一個常量。而常量對一個變數求導，高中數學告訴我們結果當然是0，所以我們返回一個新的常量`GConstant(0)`,而它的`expression`也很簡單，就是返回本身的值。

接下來是Variable

class GVariable(GBaseClass): def __init__(self, name, value=None): self.name = name self.value = value self.type = "VARIABLE" def partialGradient(self, partial): if partial.name == self.name: return GConstant(1) return GConstant(0) def expression(self): return str(self.name)

甚至比常量還簡單一些，因為是變數，所以它的值可能是不確定的，構造的時候默認為None，一個變數它對自身的導數是1，對其它變數是0，所以我們可以看到在`partialGradient()`也正是這樣操作的。變數本身的`expression`也就是它本身的標識符。

緊接著就是大頭了,`Operation`,比如圖1所示，我們將一個變數和一個常量通過二元運算元`plus`連接起來，本身它的結果就構成了一個計算式了。

class GOperation(GBaseClass): def __init__(self, a, b, operType): self.operatorType = operType; self.left = a self.right = b

幾乎所有計算都是二元的，所以我們可以傳入兩個運算元，operType是一個字元串，指示用什麼計算項連接兩個運算元。對於特殊的比如`exp`等一元計算項，可以默認傳入的右運算元為None。

接下來我們需要求偏導和寫expression了。

def partialGradient(self, partial): # partial must be a variable if partial.type != "VARIABLE": return None if self.operatorType == "plus" return GOperationWrapper(self.left.partialGradient(partial), self.right.partialGradient(partial),"plus") def expression(self): if self.operatorType == "plus": return self.left.expression() + "+" + self.right.expression()

比如我們先看看最簡單的「加法」，`GOperationWrapper`是對GOperation的外層封裝，後面一些優化可以在裡面完成，現在你可以直接認為：

def GOperationWrapper(left, right, operType): return GOperation(left, right, operType)

#求導

我們來看看`partialGradient`做了什麼，回憶一下高中數學，對一個加式的求導，就是左右兩邊運算元分別求導再相加，所以我們在`partialGradient`就翻譯了這個操作而已，複雜的事情交給遞歸去解決，`expression`同理，更加簡單。

當然此時我們只有`plus`這一個計算項，肯定無法處理複雜的情況，所以我們添加更多的計算項就可以了：

def partialGradient(self, partial): # partial must be a variable if partial.type != "VARIABLE": return None if self.operatorType == "plus" or self.operatorType == "minus": return GOperationWrapper(self.left.partialGradient(partial), self.right.partialGradient(partial), self.operatorType) if self.operatorType == "multiple": part1 = GOperationWrapper(self.left.partialGradient(partial), self.right, "multiple") part2 = GOperationWrapper(self.left, self.right.partialGradient(partial), "multiple") return GOperationWrapper(part1, part2, "plus") if self.operatorType == "division": part1 = GOperationWrapper(self.left.partialGradient(partial), self.right, "multiple") part2 = GOperationWrapper(self.left, self.right.partialGradient(partial), "multiple") part3 = GOperationWrapper(part1, part2, "minus") part4 = GOperationWrapper(self.right, GConstant(2), pow) part5 = GOperationWrapper(part3, part4, division) return part5 # pow should be g^a,a is a constant. if self.operatorType == "pow": c = GConstant(self.right.value - 1) part2 = GOperationWrapper(self.left, c, "pow") part3 = GOperationWrapper(self.right, part2, "multiple") return GOperationWrapper(self.left.partialGradient(partial), part3, "multiple") if self.operatorType == "exp": return GOperationWrapper(self.left.partialGradient(partial), self, "multiple") if self.operatorType == "ln": part1 = GOperationWrapper(GConstant(1),self.left,"division") rst = GOperationWrapper(self.left.partialGradient(partial), part1, "multiple") return rst return None

咱一個一個看：`minus`和`plus`類似，你也可以把高中課本的求導公式翻出來一個一個對照：

當然還有最最重要的鏈式法則

比如就拿稍顯複雜的`division`計算項來說：

if self.operatorType == "division": part1 = GOperationWrapper(self.left.partialGradient(partial), self.right, "multiple") part2 = GOperationWrapper(self.left, self.right.partialGradient(partial), "multiple") part3 = GOperationWrapper(part1, part2, "minus") part4 = GOperationWrapper(self.right, GConstant(2), pow) part5 = GOperationWrapper(part3, part4, division) return part5

對應的求導公式是

代碼里的`part1`就是uv,`part2`是uv,`part3`是uv-uv，`part4`是v^2,最後的結果`part5`則是除法計算項將uv-uv和v^2連接起來。**代碼做的不過是如實翻譯公式**而已。

另一個很重要的就是鏈式法則，比如我們在對`power`計算項求導的時候，（這裡限制了指數位置必須是常數），除了翻譯公式$y=x^n,y=nx^{n-1}$外，**還要考慮底數部分可能是一個函數，所以還需要乘上這個函數的偏導**：

if self.operatorType == "pow": c = GConstant(self.right.value - 1) part2 = GOperationWrapper(self.left, c, "pow") part3 = GOperationWrapper(self.right, part2, "multiple") return GOperationWrapper(self.left.partialGradient(partial), part3, "multiple")

至此我們已經完成了主要的部分，我們可以在這些基礎計算項的基礎上**封裝出更複雜的計算邏輯**，比如神經網路中常用的Sigmoid

def sigmoid(X): a = GConstant(1.0) b = GOperationWrapper(GConstant(0), X, minus) c = GOperationWrapper(b, None, exp) d = GOperationWrapper(a, c, plus) rst = GOperationWrapper(a, d, division) return rst

你完全不用關心如果對sigmoid求導，因為你只需要對它返回的結果調用`partialGradient()`就可以了，遞歸會自動去梳理其中的拓撲序，完成導數求解。

##驗證

我們試著構建一個計算式然後運行一下（完整代碼見[代碼1][3]）：

# case 3 X = GVariable("x") y = GVariable("y") beta = GVariable("beta") xb = GOperationWrapper(X, beta, multiple) s_xb = sigmoid(xb) m = GOperationWrapper(s_xb, y, minus) f = GOperationWrapper(m, GConstant(2), pow) print "F:
", f.expression() print "F partial gradient of B:
", f.partialGradient(x).expression()

上面我們構造了如下公式

程序輸出為：

F:

((1.0)/(1.0+exp(0-(x)*(beta)))-y)^(2)

F partial gradient of x:

(((0)*(1.0+exp(0-(x)*(beta)))-(1.0)*(0+(0-(1)*(beta)+(x)*(0))*(exp(0-(x)*(beta)))))/((1.0+exp(0-(x)*(beta)))^(2))-0)*((2)*(((1.0)/(1.0+exp(0-(x)*(beta)))-y)^(1)))

天啦嚕!!(╯ - )╯︵ ┻━┻ ，怎麼是這麼複雜的一堆，如何驗證結果是對的呢，你可以把上面的式子拷貝到wolframe alpha上，第一個式子的結果里我們發現wolframe alpha已經自動幫我們對x求了一次導：

第二個求導結果放進去，發現它在「alternate form」里有一個形態稍加轉化就是上面這個求導結果（分子提取一個-2出來）：

所以我們的求導結果是對的。

接下來有個問題，我們列印出來的東西太複雜了，明細有很多地方可以簡化，比如`0*a=0`、`1*a=a`這樣的小學知識就可以幫到我們，可以明顯幫我們簡化公式，這個時候就到了我們的`GOperationWrapper`了，加入一些簡單的邏輯：

def GOperationWrapper(left, right, operType): if operType == "multiple": if left.type == "CONSTANT" and right.type == "CONSTANT": return GConstant(left.value * right.value) if left.type == "CONSTANT" and left.value == 1: return right if left.type == "CONSTANT" and left.value == 0: return GConstant(0) if right.type == "CONSTANT" and right.value == 1: return left if right.type == "CONSTANT" and right.value == 0: return GConstant(0) if operType == "plus": if left.type == "CONSTANT" and left.value == 0: return right if right.type == "CONSTANT" and right.value == 0: return left if operType == "minus": if right.type == "CONSTANT" and right.value == 0: return left return GOperation(left, right, operType)

都是小學課本如實翻譯，就可以把結果簡化掉,可以看到已經減少了一截了，而且對於計算也有一些優化。完整代碼見代碼2：

F partial gradient of x:

((0-(0-beta)*(exp(0-(x)*(beta))))/((1.0+exp(0-(x)*(beta)))^(2)))*((2)*(((1.0)/(1.0+exp(0-(x)*(beta)))-y)^(1)))

#還能做什麼，優化！

接下來我們還能做什麼呢？在寫一個類似的遞歸函數傳入Variable的值然後計算函數式的結果，這個就不在這寫了，大同小異。

我們梳理下剛才調用的邏輯，你會發現對x求導到最底層的時候做了很多**重複計算**，大家回憶一下遞歸的好處，其中有一個就是「記憶化搜索」，可以大幅提高運行效率。也就是在第一次運行的時候記錄下結果，以後再調用的時候就直接返回存好的結果。

所以我們可以在 求導/求表達式 的時候把結果存下來：

比如對expression進行改造：代碼見[代碼][7]

def expression(self): if self.expressionRst != None: return self.expressionRst .... .... .... self.expressionRst = rst return rst

除此之外還可以做更多的優化，比如在不同地方可能會出現相同的計算式，其實完全可以根據計算式的expression，進一步記憶化，保證每一個式子只在程序里出現一次，**比如我們在過程中多次使用到了`GConstant(0)`，其實這個完全只聲明並使用一次**。

通過列印`G_STATIC_VARIABLE`我們發現程序運行一次創建了13個常量，而對`GConstant`進行一層記憶化封裝之後：

G_CONSTANT_DICT = {} def GConstantWrapper(value): global G_CONSTANT_DICT if G_CONSTANT_DICT.has_key(value): return G_CONSTANT_DICT[value] rst = GConstant(value) G_CONSTANT_DICT.setdefault(value,rst) return rst

最後一共只創建了3個常量，(0)，(1)和(2)，這些東西都可以重複利用，不需要浪費空間和CPU去聲明新實例，這也符合函數式編程的思想，這這裡推薦大家讀一下《SICP》，會有幫助的。**我們甚至可以將這個思路推廣到所有出現的計算式**，可以在後續計算和求導的時候節省大量的時間，不過在此就不做實現了。

#尾巴

花了一個小時寫代碼，N個碎片時間寫博客，但真心覺得求導鏈式法則和遞歸簡直就是天作之合，不記錄一下於心難忍。當然真實tf和mxnet使用的自動求導肯定還有更多優化的，不過就不深鑽下去了，這個狀態~味道剛剛好。