作者丨蘇劍林
單位丨廣州火焰資訊科技有限公司
研究方向丨NLP，神經網路
個人主頁丨kexue.fm
本來筆者已經決心不玩 RNN 了，但是在上個星期思考時忽然意識到 RNN 實際上對應了 ODE（常微分方程）的數值解法，這為我一直以來想做的事情——用深度學習來解決一些純數學問題——提供了思路。事實上這是一個頗為有趣和有用的結果，遂介紹一翻。順便地，本文也涉及到了自己動手編寫 RNN 的內容，所以本文也可以作為編寫自定義的 RNN 層的一個簡單教程。
注：本文並非前段時間的熱點“神經 ODE [1]”的介紹（但有一定的聯絡）。
RNN基本
什麼是RNN？
眾所周知，RNN 是“迴圈神經網路（Recurrent Neural Network）”，跟 CNN 不同，RNN 可以說是一類模型的總稱，而並非單個模型。簡單來講，只要是輸入向量序列 (x1,x2,…,xT)，輸出另外一個向量序列 (y1,y2,…,yT)，並且滿足如下遞迴關係的模型，都可以稱為 RNN。
也正因為如此，原始的樸素 RNN，還有改進的如 GRU、LSTM、SRU 等模型，我們都稱為 RNN，因為它們都可以作為上式的一個特例。還有一些看上去與 RNN 沒關的內容，比如前不久介紹的 CRF 的分母的計算，實際上也是一個簡單的 RNN。
說白了，RNN 其實就是遞迴計算。
自己編寫RNN
這裡我們先介紹如何用 Keras 簡單快捷地編寫自定義的 RNN。
事實上，不管在 Keras 還是純 tensorflow 中，要自定義自己的 RNN 都不算複雜。在 Keras 中，只要寫出每一步的遞迴函式；而在 tensorflow 中，則稍微複雜一點，需要將每一步的遞迴函式封裝為一個 RNNCell 類。
下面介紹用 Keras 實現最基本的一個 RNN：
程式碼非常簡單：
#! -*- coding: utf-8- -*-
from keras.layers import Layer
import keras.backend as K
class My_RNN(Layer):
def __init__(self, output_dim, **kwargs):
self.output_dim = output_dim # 輸出維度
super(My_RNN, self).__init__(**kwargs)
def build(self, input_shape): # 定義可訓練引數
self.kernel1 = self.add_weight(name='kernel1',
shape=(self.output_dim, self.output_dim),
initializer='glorot_normal',
trainable=True)
self.kernel2 = self.add_weight(name='kernel2',
shape=(input_shape[-1], self.output_dim),
initializer='glorot_normal',
trainable=True)
self.bias = self.add_weight(name='kernel',
shape=(self.output_dim,),
initializer='glorot_normal',
trainable=True)
def step_do(self, step_in, states): # 定義每一步的迭代
step_out = K.tanh(K.dot(states[0], self.kernel1) +
K.dot(step_in, self.kernel2) +
self.bias)
return step_out, [step_out]
def call(self, inputs): # 定義正式執行的函式
init_states = [K.zeros((K.shape(inputs)[0],
self.output_dim)
)] # 定義初始態(全零)
outputs = K.rnn(self.step_do, inputs, init_states) # 迴圈執行step_do函式
return outputs[0] # outputs是一個tuple，outputs[0]為最後時刻的輸出，
# outputs[1]為整個輸出的時間序列，output[2]是一個list，
# 是中間的隱藏狀態。
def compute_output_shape(self, input_shape):
return (input_shape[0], self.output_dim)可以看到，雖然程式碼行數不少，但大部分都只是固定格式的語句，真正定義 RNN 的，是 step_do 這個函式，這個函式接受兩個輸入：step_in 和 states。其中 step_in 是一個 (batch_size, input_dim) 的張量，代表當前時刻的樣本 xt，而 states 是一個 list，代表 yt−1 及一些中間變數。
特別要注意的是，states 是一個張量的 list，而不是單個張量，這是因為在遞迴過程中可能要同時傳遞多箇中間變數，而不僅僅是 yt−1 一個，比如 LSTM 就需要有兩個態張量。最後 step_do 要返回 yt 和新的 states，這是 step_do 這步的函式的編寫規範。
而 K.rnn 這個函式，接受三個基本引數（還有其他引數，請自行看官方文件），其中第一個引數就是剛才寫好的 step_do 函式，第二個引數則是輸入的時間序列，第三個是初始態，跟前面說的 states 一致，所以很自然 init_states 也是一個張量的 list，預設情況下我們會選擇全零初始化。
ODE基本
什麼是ODE？
ODE 就是“常微分方程（Ordinary Differential Equation）”，這裡指的是一般的常微分方程組：
研究 ODE 的領域通常也直接稱為“動力學”、“動力系統”，這是因為牛頓力學通常也就只是一組 ODE 而已。
ODE可以產生非常豐富的函式。比如 e^t 其實就是 x˙=x 的解，sint 和 cost 都是 x¨+x=0 的解（初始條件不同）。事實上，我記得確實有一些教程是直接通過微分方程 x˙=x 來定義 e^t 函式的。除了這些初等函式，很多我們能叫得上名字但不知道是什麼鬼的特殊函式，都是通過 ODE 匯出來的，比如超幾何函式、勒讓德函式、貝塞爾函式...
總之，ODE 能產生並且已經產生了各種各樣千奇百怪的函式～
數值解ODE
能精確求出解析解的 ODE 其實是非常少的，所以很多時候我們都需要數值解法。
ODE 的數值解已經是一門非常成熟的學科了，這裡我們也不多做介紹，僅引入最基本的由數學家尤拉提出來的迭代公式：
這裡的 h 是步長。尤拉的解法來源很簡單，就是用：
來近似導數項 x˙(t)。只要給定初始條件 x(0)，我們就可以根據 (4) 一步步迭代算出每個時間點的結果。
ODE與RNN
ODE也是RNN
大家仔細對比 (4) 和 (1)，發現有什麼聯絡了嗎？
在 (1) 中，t 是一個整數變數，在 (4) 中，t 是一個浮點變數，除此之外，(4) 跟 (1) 貌似就沒有什麼明顯的區別了。事實上，在 (4) 中我們可以以 h 為時間單位，記 t=nh，那麼 (4) 變成了：
可以看到現在 (6) 中的時間變數 n 也是整數了。這樣一來，我們就知道：ODE 的尤拉解法 (4) 實際上就是 RNN 的一個特例罷了。這裡我們也許可以間接明白為什麼 RNN 的擬合能力如此之強了（尤其是對於時間序列資料），我們看到 ODE 可以產生很多複雜的函式，而 ODE 只不過是 RNN 的一個特例罷了，所以 RNN 也就可以產生更為複雜的函式了。
用RNN解ODE
於是，我們就可以寫一個 RNN 來解 ODE 了，比如《兩生物種群競爭模型》[2] 中的例子：
我們可以寫出：
#! -*- coding: utf-8- -*-
from keras.layers import Layer
import keras.backend as K
class ODE_RNN(Layer):
def __init__(self, steps, h, **kwargs):
self.steps = steps
self.h = h
super(ODE_RNN, self).__init__(**kwargs)
def step_do(self, step_in, states): # 定義每一步的迭代
x = states[0]
r1,r2,a1,a2,iN1,iN2 = 0.1,0.3,0.0001,0.0002,0.002,0.003
_1 = r1 * x[:,0] * (1 - iN1 * x[:,0]) - a1 * x[:,0] * x[:,1]
_2 = r2 * x[:,1] * (1 - iN2 * x[:,1]) - a2 * x[:,0] * x[:,1]
_1 = K.expand_dims(_1, 1)
_2 = K.expand_dims(_2, 1)
_ = K.concatenate([_1, _2], 1)
step_out = x + self.h * _
return step_out, [step_out]
def call(self, inputs): # 這裡的inputs就是初始條件
init_states = [inputs]
zeros = K.zeros((K.shape(inputs)[0],
self.steps,
K.shape(inputs)[1])) # 迭代過程用不著外部輸入，所以
# 指定一個全零輸入，只為形式上的傳入
outputs = K.rnn(self.step_do, zeros, init_states) # 迴圈執行step_do函式
return outputs[1] # 這次我們輸出整個結果序列
def compute_output_shape(self, input_shape):
return (input_shape[0], self.steps, input_shape[1])
from keras.models import Sequential
import numpy as np
import matplotlib.pyplot as plt
steps,h = 1000,0.1
M = Sequential()
M.add(ODE_RNN(steps, h, input_shape=(2,)))
M.summary()
# 直接前向傳播就輸出解了
result = M.predict(np.array([[100, 150]]))[0] # 以[100, 150]為初始條件進行演算
times = np.arange(1, steps+1) * h
# 繪圖
plt.plot(times, result[:,0])
plt.plot(times, result[:,1])
plt.savefig('test.png')整個過程很容易理解，只不過有兩點需要指出一下。首先，由於方程組 (7) 只有兩維，而且不容易寫成矩陣運算，因此我在 step_do 函式中是直接逐位操作的（程式碼中的 x[:,0],x[:,1]），如果方程本身維度較高，而且能寫成矩陣運算，那麼直接利用矩陣運算寫會更加高效；然後，我們可以看到，寫完整個模型之後，直接 predict 就輸出結果了，不需要“訓練”。
▲ RNN解兩物種的競爭模型
反推ODE引數
前一節的介紹也就是說，RNN 的前向傳播跟 ODE 的尤拉解法是對應的，那麼反向傳播又對應什麼呢？
在實際問題中，有一類問題稱為“模型推斷”，它是在已知實驗資料的基礎上，猜測這批資料符合的模型（機理推斷）。這類問題的做法大概分兩步，第一步是猜測模型的形式，第二步是確定模型的引數。假定這批資料可以由一個 ODE 描述，並且這個 ODE 的形式已經知道了，那麼就需要估計裡邊的引數。
如果能夠用公式完全解出這個 ODE，那麼這就只是一個非常簡單的迴歸問題罷了。但前面已經說過，多數 ODE 都沒有公式解，所以數值方法就必須了。這其實就是 ODE 對應的 RNN 的反向傳播所要做的事情：前向傳播就是解 ODE（RNN 的預測過程），反向傳播自然就是推斷 ODE 的引數了（RNN 的訓練過程）。這是一個非常有趣的事實：ODE 的引數推斷是一個被研究得很充分的課題，然而在深度學習這裡，只是 RNN 的一個最基本的應用罷了。
我們把剛才的例子的微分方程的解資料儲存下來，然後只取幾個點，看看能不能反推原來的微分方程出來，解資料為：
假設就已知這有限的點資料，然後假定方程 (7) 的形式，求方程的各個引數。我們修改一下前面的程式碼：
#! -*- coding: utf-8- -*-
from keras.layers import Layer
import keras.backend as K
def my_init(shape, dtype=None): # 需要定義好初始化，這相當於需要實驗估計引數的量級
return K.variable([0.1, 0.1, 0.001, 0.001, 0.001, 0.001])
class ODE_RNN(Layer):
def __init__(self, steps, h, **kwargs):
self.steps = steps
self.h = h
super(ODE_RNN, self).__init__(**kwargs)
def build(self, input_shape): # 將原來的引數設為可訓練的引數
self.kernel = self.add_weight(name='kernel',
shape=(6,),
initializer=my_init,
trainable=True)
def step_do(self, step_in, states): # 定義每一步的迭代
x = states[0]
r1,r2,a1,a2,iN1,iN2 = (self.kernel[0], self.kernel[1],
self.kernel[2], self.kernel[3],
self.kernel[4], self.kernel[5])
_1 = r1 * x[:,0] * (1 - iN1 * x[:,0]) - a1 * x[:,0] * x[:,1]
_2 = r2 * x[:,1] * (1 - iN2 * x[:,1]) - a2 * x[:,0] * x[:,1]
_1 = K.expand_dims(_1, 1)
_2 = K.expand_dims(_2, 1)
_ = K.concatenate([_1, _2], 1)
step_out = x + self.h * K.clip(_, -1e5, 1e5) # 防止梯度爆炸
return step_out, [step_out]
def call(self, inputs): # 這裡的inputs就是初始條件
init_states = [inputs]
zeros = K.zeros((K.shape(inputs)[0],
self.steps,
K.shape(inputs)[1])) # 迭代過程用不著外部輸入，所以
# 指定一個全零輸入，只為形式上的傳入
outputs = K.rnn(self.step_do, zeros, init_states) # 迴圈執行step_do函式
return outputs[1] # 這次我們輸出整個結果序列
def compute_output_shape(self, input_shape):
return (input_shape[0], self.steps, input_shape[1])
from keras.models import Sequential
from keras.optimizers import Adam
import numpy as np
import matplotlib.pyplot as plt
steps,h = 50, 1 # 用大步長，減少步數，削弱長時依賴，也加快推斷速度
series = {0: [100, 150],
10: [165, 283],
15: [197, 290],
30: [280, 276],
36: [305, 269],
40: [318, 266],
42: [324, 264]}
M = Sequential()
M.add(ODE_RNN(steps, h, input_shape=(2,)))
M.summary()
# 構建訓練樣本
# 其實就只有一個樣本序列，X為初始條件，Y為後續時間序列
X = np.array([series[0]])
Y = np.zeros((1, steps, 2))
for i,j in series.items():
if i != 0:
Y[0, int(i/h)-1] += series[i]
# 自定義loss
# 在訓練的時候，只考慮有資料的幾個時刻，沒有資料的時刻被忽略
def ode_loss(y_true, y_pred):
T = K.sum(K.abs(y_true), 2, keepdims=True)
T = K.cast(K.greater(T, 1e-3), 'float32')
return K.sum(T * K.square(y_true - y_pred), [1, 2])
M.compile(loss=ode_loss,
optimizer=Adam(1e-4))
M.fit(X, Y, epochs=10000) # 用低學習率訓練足夠多輪
# 用訓練出來的模型重新預測，繪圖，比較結果
result = M.predict(np.array([[100, 150]]))[0]
times = np.arange(1, steps+1) * h
plt.clf()
plt.plot(times, result[:,0], color='blue')
plt.plot(times, result[:,1], color='green')
plt.plot(series.keys(), [i[0] for i in series.values()], 'o', color='blue')
plt.plot(series.keys(), [i[1] for i in series.values()], 'o', color='green')
plt.savefig('test.png')結果可以用一張圖來看：
▲ RNN做ODE的引數估計效果
（散點：有限的實驗資料，曲線：估計出來的模型）
顯然結果是讓人滿意的。
又到總結
本文在一個一般的框架下介紹了 RNN 模型及其在 Keras 下的自定義寫法，然後揭示了 ODE 與 RNN 的聯絡。在此基礎上，介紹了用 RNN 直接求解 ODE 以及用 RNN 反推 ODE 引數的基本思路。
需要提醒讀者的是，在 RNN 模型的反向傳播中，要謹慎地做好初始化和截斷處理處理，並且選擇好學習率等，以防止梯度爆炸的出現（梯度消失只是優化得不夠好，梯度爆炸則是直接崩潰了，解決梯度爆炸問題尤為重要）。
總之，梯度消失和梯度爆炸在 RNN 中是一個很經典的困難，事實上，LSTM、GRU 等模型的引入，根本原因就是為了解決 RNN 的梯度消失問題，而梯度爆炸則是通過使用 tanh 或 sigmoid 啟用函式來解決的。
但是如果用 RNN 解決 ODE 的話，我們就沒有選擇啟用函式的權利了（啟用函式就是 ODE 的一部分），所以只能謹慎地做好初始化及其他處理。據說，只要謹慎做好初始化，普通 RNN 中用 relu 作為啟用函式都是無妨的。
相關連結
[1]. Tian Qi C, Yulia R, Jesse B, David D. Neural Ordinary Differential Equations. arXiv preprint arXiv:1806.07366, 2018.
[2]. 兩生物種群競爭模型
https://kexue.fm/archives/3120