【Python機器學習】用遺傳演算法實現符號回歸——淺析gplearn

01-27

整理自：Introduction to GP

提起Python里的機器學習包，我們通常會想到Scikit-Learn、XGBoost、LightGBM、Keras和TensorFlow。在這裡，我將向大家介紹gplearn。它是目前Python內最成熟的符號回歸演算法實現。

簡介

在線性模型中，目標變數y可以表示為 $y=textbf{w}^Ttextbf{x}+b$ 。然而，即使加入了多項式特徵，它也很難發現特徵變數與目標變數之間更複雜的關係。

作為一種一種監督學習方法，符號回歸（symbolic regression）試圖發現某種隱藏的數學公式，以此利用特徵變數預測目標變數。

符號回歸的具體實現方式是遺傳演算法（genetic algorithm）。一開始，一群天真而未經歷選擇的公式會被隨機生成。此後的每一代中，最「合適」的公式們將被選中。隨著迭代次數的增長，它們不斷繁殖、變異、進化，從而不斷逼近數據分布的真相。

公式的表示方法

假設我們有特徵X0和X1，需要預測目標y。一個可能的公式是：

$y = X_0^{2} - 3 times X_1 + 0.5$

它也可以寫作：

$y = X_0 times X_0 - 3 times X_1 + 0.5$

同時，我們也可以把它轉化成一個S-表達式：

$y = (+ (- (times X_0 X_0) (times 3 X_1)) 0.5)$

公式里包括了變數（X0和X1）、函數（加、減、乘）和常數（3和0.5）。

有了S-表達式，我們可以把公式表示為一個二叉樹：

在這個二叉樹里，所有的葉節點都是變數或者常數，內部的節點則是函數；公式的輸出值可以用遞歸的方法求得。

在用遺傳演算法實現符號回歸時，所有的公式都會以S-表達式和二叉樹的形式存在。值得注意的是，上圖樹內的任意子樹都可以被修改或替代。例如，如果把 $y = (+ (- (times X_0 X_0) (times 3 X_1)) 0.5)$ 中的 $(times 3 X_1)$ 替換為 $(times 6 X_0)$ ，那麼我們就得到了一個新的公式： $y = (+ (- (times X_0 X_0) (times 6 X_0)) 0.5)$ 。

gplearn包括了一系列不同的函數，以供用戶選擇。它們的關鍵詞是：

add：加法，二元運算
sub：減法，二元運算
mul：乘法，二元運算
div：除法，二元運算
sqrt：平方根，一元運算
log：對數，一元運算
abs：絕對值，一元運算
neg：相反數，一元運算
inv：倒數，一元運算
max：最大值，二元運算
min：最小值，二元運算
sin：正弦（弧度），一元運算
cos：餘弦（弧度），一元運算
tan：正切（弧度），一元運算

因為輸入值很可能是隨機產生的，所以這些函數需要處理諸如「除以零」的問題。因此，gplearn中的許多函數並不符合它們原本的數學定義，而是「受保護」的修改版：

除法：如果 $-0.001 le b le 0.001$ ，那麼定義 ${div}_{safe} (a, b)=1$ 。
平方根：定義 ${sqrt}_{safe}(a) = sqrt {|a|}$
對數：如果 $-0.001 le a le 0.001$ ，那麼 ${log}_{safe} (a) = 0$ ，其他情況下定義 ${log}_{safe} (a) = log |a|$
倒數：如果 $-0.001 le a le 0.001$ ，那麼定義 ${inv}_{safe}(a)=1$ 。

gplearn還允許用戶自定義函數。

公式的適應度

物競天擇，適者生存。——《天演論》

和其他機器學習演算法一樣，遺傳演算法的核心在於衡量公式的適應度（fitness function）。在符號回歸里，適應度的地位類似於目標函數、score、loss和error。

gplearn的主要組成部分有兩個：SymbolicRegressor和SymbolicTransformer。兩者的適應度有所不同。

SymbolicRegressor是回歸器。它利用遺傳演算法得到的公式，直接預測目標變數的值。
SymbolicTransformer是轉換器。它並不直接預測目標變數，而是轉化原有的特徵、輸出新的特徵，這在特徵工程的階段尤為有效。

SymbolicRegressor的適應度有三種，都是機器學習里常見的error function：

mae: mean absolute error
mse: mean squared error
rmse: root mean squared error

SymbolicTransformer會最大化輸出的新特徵與目標變數之間的相關係數的絕對值：（並非相關係數本身，因為很大的負相關反而有利於預測）

pearson：皮爾遜積矩相關係數（Pearson product-moment correlation coefficient），默認設置。因為它衡量線性相關度，所以適用於下一個預測器（分類或回歸）是線性模型的情況。
spearman：斯皮爾曼等級相關係數（Spearmans rank correlation coefficient）。因為它衡量單調相關度，所以適用於下一個預測器（分類或回歸）是決策樹類模型的情況。

當然，用戶也可以自定義適應度的標準。

遺傳演算法內，耗時最大的部分無疑是適應度的計算。所以，gplearn允許用戶通過修改n_jobs參數控制並行運算。在數據量和公式數量較大時，並行計算的速度優勢最為明顯。

公式的進化

gplearn中最重要的超參數都與公式的進化方式有關。

在初始化階段，population_size棵公式樹（參見上文「公式的表達方法」）會被隨機生成。每棵公式樹的深度都會受到init_depth參數的限制。init_depth是一個二元組(min_depth, max_depth)，樹的初始深度將處在[min_depth, max_depth]的區間內（包含端點）。

通常而言，變數越多，模型越複雜，那麼population_size就越大越好。

每棵公式樹的初始化方式由init_method控制，分為三種策略：

grow：公式樹從根節點開始生長。在每一個子節點，gplearn會從所有常數、變數和函數中隨機選取一個元素。如果它是常數或者變數，那麼這個節點會停止生長，成為一個葉節點。如果它是函數，那麼它的兩個子節點將繼續生長。用grow策略生長得到的公式樹往往不對稱，而且普遍會比用戶設置的最大深度淺一些；在變數的數量遠大於函數的數量時，這種情況更明顯。
full：除了最後一層外，其他所有層的所有節點都是內部節點——它們都只是隨機選擇的函數，而不是變數或者常數。最後一層的葉節點則是隨機選擇的變數和常數。用full策略得到的公式樹必然是perfect binary tree。
half and half：一半的公式樹用grow策略生成，另一半用full策略生成。因為種群的多樣性有利於生存，所以這是init_method參數的默認值。

為了模擬自然選擇的過程，大部分「不適應環境」，即適應度不足的公式會被淘汰。從每一代的所有公式中，tournament_size個公式會被隨機選中，其中適應度最高的公式將被認定為生存競爭的勝利者，進入下一代。tournament_size的大小與進化論中的選擇壓力息息相關：tournament_size越小，選擇壓力越大，演算法收斂的速度可能更快，但也有可能錯過一些隱藏的優秀公式。

進入下一代的優勝者未必原封不動——完全不改變優勝者，直接讓它進入下一代的策略被稱為繁殖（reproduction）。用戶可以採取一系列的變異措施：

交叉（Crossover）

優勝者內隨機選擇一個子樹，替換為另一棵公式樹的隨機子樹。此處的另一棵公式樹通常是剩餘公式樹中適應度最高的。

子樹變異（Subtree Mutation）

由p_subtree_mutation參數控制。這是一種更激進的變異策略：優勝者的一棵子樹將被另一棵完全隨機的全新子樹代替。

hoist變異（Hoist Mutation）

由p_hoist_mutation參數控制。hoist變異是一種對抗公式樹膨脹（bloating，即過於複雜）的方法：從優勝者公式樹內隨機選擇一個子樹A，再從A里隨機選擇一個子樹B，然後把B提升到A原來的位置，用B替代A。hoist的含義即「升高、提起」。

點變異（Point Mutation）

由p_point_replace參數控制。一個隨機的節點將會被改變，比如加法可以被替換成除法，變數X0可以被替換成常數-2.5。點變異可以重新加入一些先前被淘汰的函數和變數，從而促進公式的多樣性。

由於進化過程的黑箱屬性，調參很大程度上依賴於用戶的直覺和經驗。對於實際問題本身的理解也必不可少，如：最終的公式可能有多複雜？

在使用SymbolicRegressor時，對目標變數進行標準化（standardization）和區間縮放（min-max-scaling）可以有效避免常數值區間不符的問題。（假如目標變數的區間是[-1000, 4000]，那麼區間為[-1, 1]的常數恐怕不會有多大幫助，最終得出的公式只會是一串幾乎毫無意義的加法和乘法。）

公式樹的膨脹以及解決辦法

一棵公式樹的複雜度有兩個方面：深度（樹的深度）和長度（節點的總數量）。當公式變得越來越複雜，計算速度也越發緩慢，但它的適應度卻毫無提升時，我們稱這種現象為膨脹（bloating）。

對抗膨脹的方法如下：

在適應度函數中加入節儉係數（parsimony coefficient），由參數parsimony_coefficient控制，懲罰過於複雜的公式。節儉係數往往由實踐驗證決定。如果過於吝嗇（節儉係數太大），那麼所有的公式樹都會縮小到只剩一個變數或常數；如果過於慷慨（節儉係數太小），公式樹將嚴重膨脹。不過，gplearn已經提供了auto的選項，能自動控制節儉項的大小。
使用hoist變異，削減過大的公式樹。
每一代的進化中，只有一部分樣本會被隨機選取，用于衡量公式的適應度。此處樣本的數量由參數max_samples控制。這種做法不僅提升了運算速度、保證了公式的多樣性，還允許用戶查看每一個公式的out-of-bag fitness，更為客觀。

例：符號回歸 vs. 決策樹 vs. 隨機森林

首先，創造一個基於 $y = X_0^{2} - X_1^{2} + X_1 - 1$ 的數據生成分布：

x0 = np.arange(-1, 1, 1/10.)nx1 = np.arange(-1, 1, 1/10.)nx0, x1 = np.meshgrid(x0, x1)ny_truth = x0**2 - x1**2 + x1 - 1nnax = plt.figure().gca(projection=3d)nax.set_xlim(-1, 1)nax.set_ylim(-1, 1)nsurf = ax.plot_surface(x0, x1, y_truth, rstride=1, cstride=1,n color=green, alpha=0.5)nplt.show()n

得出的圖象為：

然後，根據數據生成分布生成隨機的訓練集和測試集：

rng = check_random_state(0)nn# Training samplesnX_train = rng.uniform(-1, 1, 100).reshape(50, 2)ny_train = X_train[:, 0]**2 - X_train[:, 1]**2 + X_train[:, 1] - 1nn# Testing samplesnX_test = rng.uniform(-1, 1, 100).reshape(50, 2)ny_test = X_test[:, 0]**2 - X_test[:, 1]**2 + X_test[:, 1] - 1n

接下來，我們訓練一個SymbolicRegressor：

est_gp = SymbolicRegressor(population_size=5000,n generations=20, stopping_criteria=0.01,n p_crossover=0.7, p_subtree_mutation=0.1,n p_hoist_mutation=0.05, p_point_mutation=0.1,n max_samples=0.9, verbose=1,n parsimony_coefficient=0.01, random_state=0)nest_gp.fit(X_train, y_train)nn | Population Average | Best Individual |n---- ------------------------- ------------------------------------------ ----------n Gen Length Fitness Length Fitness OOB Fitness Time Leftn 0 38.13 386.19117972 7 0.331580808730 0.470286152255 55.15sn 1 9.91 1.66832489614 5 0.335361761359 0.488347149514 1.25mn 2 7.76 1.888657267 7 0.260765934398 0.565517599814 1.45mn 3 5.37 1.00018638338 17 0.223753461954 0.274920433701 1.42mn 4 4.69 0.878161643513 17 0.145095322600 0.158359554221 1.35mn 5 6.1 0.91987274474 11 0.043612562970 0.043612562970 1.31mn 6 7.18 1.09868887802 11 0.043612562970 0.043612562970 1.23mn 7 7.65 1.96650325011 11 0.043612562970 0.043612562970 1.18mn 8 8.02 1.02643443398 11 0.043612562970 0.043612562970 1.08mn 9 9.07 1.22732144371 11 0.000781474035 0.0007814740353 59.43sn

作為比較，我們同時使用scikit-learn的決策樹和隨機森林進行訓練：

est_tree = DecisionTreeRegressor()nest_tree.fit(X_train, y_train)nest_rf = RandomForestRegressor()nest_rf.fit(X_train, y_train)nny_gp = est_gp.predict(np.c_[x0.ravel(), x1.ravel()]).reshape(x0.shape)nscore_gp = est_gp.score(X_test, y_test)ny_tree = est_tree.predict(np.c_[x0.ravel(), x1.ravel()]).reshape(x0.shape)nscore_tree = est_tree.score(X_test, y_test)ny_rf = est_rf.predict(np.c_[x0.ravel(), x1.ravel()]).reshape(x0.shape)nscore_rf = est_rf.score(X_test, y_test)nnfig = plt.figure(figsize=(12, 10))nnfor i, (y, score, title) in enumerate([(y_truth, None, "Ground Truth"),n (y_gp, score_gp, "SymbolicRegressor"),n (y_tree, score_tree, "DecisionTreeRegressor"),n (y_rf, score_rf, "RandomForestRegressor")]):nn ax = fig.add_subplot(2, 2, i+1, projection=3d)n ax.set_xlim(-1, 1)n ax.set_ylim(-1, 1)n surf = ax.plot_surface(x0, x1, y, rstride=1, cstride=1, color=green, alpha=0.5)n points = ax.scatter(X_train[:, 0], X_train[:, 1], y_train)n if score is not None:n score = ax.text(-.7, 1, .2, "$R^2 =/ %.6f$" % score, x, fontsize=14)n plt.title(title)nplt.show()n

作圖的結果如下：

如圖所示，SymbolicRegressor幾乎完美地還原了數據生成分布。相比之下，決策樹和隨機森林的決策邊界粗糙了不少。由此可見，已知數據生成分布很有可能可以表示為數學公式的情況下，使用符號回歸是一個優秀的選擇。

輸出公式

我們並不希望每次使用同一個公式時都要重新訓練一遍模型。恰恰相反，符號回歸的強大之處在於，我們可以保存每次訓練得到的公式。

SymbolicTransformer和SymbolicRegressor都重載了python的print函數。列印一個SymbolicRegressor時，我們會得到最優的回歸公式。列印一個SymbolicTransformer時，我們會得到一個包括了n_components個最優轉換公式的列表。

為了方便演示，我們用一個隨機生成的數據集進行訓練。

import numpy as npnfrom gplearn import geneticnnm1 = genetic.SymbolicTransformer(verbose=1, generations=3)nm1.fit(np.random.rand(10,5), np.random.rand(10))nnm2 = genetic.SymbolicRegressor(verbose=1, generations=3)nm2.fit(np.random.rand(10,5), np.random.rand(10))n

列印的結果為：

>>> print(m1)n[mul(0.180, 0.908),n mul(div(add(X4, X1), sub(X0, sub(X1, X0))), X2),n div(add(X1, X4), add(sub(X4, div(X1, sub(X1, X0))), sub(add(mul(X3, X3), X2), X4))),n div(sub(X2, div(X3, 0.278)), sub(X2, X0)),n mul(div(sub(sub(X1, 0.481), mul(0.180, 0.908)), sub(X0, X1)), X2),n div(X1, sub(X3, X1)),n add(div(mul(mul(mul(X1, X1), add(X3, X1)), div(add(X1, X3), sub(X3, -0.058))), add(mul(add(X1, X0), sub(X4, -0.003)), sub(add(X4, X1), mul(-0.273, -0.567)))), add(div(add(add(X1, X4), mul(X1, X3)), add(mul(X1, X3), sub(X2, X0))), sub(add(add(X0, X1), div(0.575, X3)), add(add(X2, X3), sub(X2, X0))))),n div(mul(div(-0.347, X1), div(X0, X3)), div(mul(0.180, 0.908), sub(X2, X0))),n div(sub(X3, X0), sub(X1, sub(add(mul(X3, 0.713), div(X1, sub(X1, X0))), mul(mul(0.164, X2), div(X3, X0))))),n add(div(div(0.575, X3), add(mul(add(X1, X0), sub(X4, -0.003)), sub(add(-0.120, X0), mul(-0.273, -0.567)))), add(div(add(add(X1, X4), mul(X1, X3)), add(mul(X1, X3), sub(X3, X1))), sub(add(add(X0, X1), div(0.575, X3)), add(add(X2, X3), sub(X2, X0)))))]n>>> print(m2)nadd(0.452, mul(mul(X1, X2), sub(X1, X3)))n

其中，X0, X2, X3等變數分別對應numpy array中的X[:, 0]、X[:, 2]和X[:, 3]。

安裝

gplearn的安裝過程很簡單，只要已有scikit-learn、numpy和scipy等常用工具包，在命令行輸入以下指令即可：

pip install gplearnn

如果需要安裝最新的開發版，可以輸入：

git clone https://github.com/trevorstephens/gplearn.gitn

然後前往gplearn所在文件夾，輸入：

python setup.py installn