内容

• ライブラリはただの便利ツールでブラックボックスなのか
• 開発におけるツール選定の本質
• 機械学習専用ライブラリと自作モデルではどこまで性能の違いが出るのか

これらを、自作したモデルと機械学習ライブラリsklearn、lightGBMを使って比較を行いました。

▼ 検証用のGBDTモデル

from sklearn.ensemble import GradientBoostingRegressor
import time

def train_sklearn_gbdt(X, y, n_trees=10, lr=0.1):
    model = GradientBoostingRegressor(
        n_estimators=n_trees,
        learning_rate=lr,
        max_depth=1,        # 自作コードと合わせる（深さ1）
        subsample=1.0,
        criterion="squared_error"
    )

    start = time.perf_counter()
    model.fit(X, y)
    end = time.perf_counter()

    print("sklearn GBDT 学習時間:", (end - start)/60)
    return model

import lightgbm as lgb
import time

def train_lightgbm(X, y, n_trees=10, lr=0.1):
    train_data = lgb.Dataset(X, label=y)

    params = {
        "objective": "regression",
        "learning_rate": lr,
        "num_leaves": 2,     # 深さ1の木に相当
        "max_depth": 1,
        "min_data_in_leaf": 1,
        "verbose": -1
    }

    start = time.perf_counter()
    model = lgb.train(params, train_data, num_boost_round=n_trees)
    end = time.perf_counter()

    print("LightGBM 学習時間:", (end - start)/60)
    return model

import numpy as np
import time

def compute_grad_hess(y, y_pred):
    grad = y_pred - y          # 勾配
    hess = np.ones_like(y)     # 2乗誤差では常に1
    return grad, hess

def calc_gain(G_L, H_L, G_R, H_R, lam=1.0):
    gain = (G_L**2 / (H_L + lam)) + (G_R**2 / (H_R + lam)) \
           - ((G_L + G_R)**2 / (H_L + H_R + lam))
    return gain

def find_best_split(X, grad, hess):
    n_samples, n_features = X.shape
    best_gain = -np.inf
    best_feature = None
    best_threshold = None

    for f in range(n_features):
        values = X[:, f]
        sorted_idx = np.argsort(values)
        values_sorted = values[sorted_idx]
        grad_sorted = grad[sorted_idx]
        hess_sorted = hess[sorted_idx]

        G_total = grad_sorted.sum()
        H_total = hess_sorted.sum()

        G_L = 0.0
        H_L = 0.0

        for i in range(1, n_samples):
            G_L += grad_sorted[i-1]
            H_L += hess_sorted[i-1]

            G_R = G_total - G_L
            H_R = H_total - H_L

            if values_sorted[i] == values_sorted[i-1]:
                continue

            gain = calc_gain(G_L, H_L, G_R, H_R)

            if gain > best_gain:
                best_gain = gain
                best_feature = f
                best_threshold = (values_sorted[i] + values_sorted[i-1]) / 2

    return best_feature, best_threshold, best_gain

def compute_leaf_value(grad, hess, lam=1.0):
    return -grad.sum() / (hess.sum() + lam)

def train_one_tree(X, y, y_pred):
    grad, hess = compute_grad_hess(y, y_pred)
    f, thr, gain = find_best_split(X, grad, hess)

    left_idx = X[:, f] <= thr
    right_idx = ~left_idx

    leaf_left = compute_leaf_value(grad[left_idx], hess[left_idx])
    leaf_right = compute_leaf_value(grad[right_idx], hess[right_idx])

    return f, thr, leaf_left, leaf_right

def train_gbdt(X, y, n_trees=10, lr=0.1):

    n_samples = len(y)
    y_pred = np.full(n_samples, y.mean()) 

    trees = []
    start = time.perf_counter()
    for t in range(n_trees):
        f, thr, leaf_L, leaf_R = train_one_tree(X, y, y_pred)

        left_idx = X[:, f] <= thr
        right_idx = ~left_idx

        y_pred[left_idx] += lr * leaf_L
        y_pred[right_idx] += lr * leaf_R

        trees.append((f, thr, leaf_L, leaf_R))
    end = time.perf_counter()
    print("GBDT 学習時間:", (end - start)/60)
    return trees, y_pred

※自作モデルに関する補足※

今回使用するモデルはGBDTの回帰問題ということで、勾配・gain・分岐・予測値更新の計算式は前回紹介したものと同じです。

検証

学習に使うデータは0～1万までの乱数を用意し、適度な計算負荷を与えました。また条件を統一するために、正則化項と学習率は一般的な値を使用し、パラメーターを調整したうえで学習を行いました。

from sklearn.metrics import mean_squared_error, mean_absolute_error
import numpo as np

# 学習データの用意
data=np.random.uniform(0,10000,(10000,3))
X=data[:,:2]
y=data[:,2]

# 自作GBDT
trees, pred = train_gbdt(X, y, n_trees=10)

# sklearn
model_sklearn = train_sklearn_gbdt(X, y, n_trees=10)
y_pred_sklearn = model_sklearn.predict(X)

# LightGBM
model_lgb = train_lightgbm(X, y, n_trees=10)
y_pred_lgb = model_lgb.predict(X)


print("GBDT MSE:", mean_squared_error(y, pred))
print("sklearn MSE:", mean_squared_error(y, y_pred_sklearn))
print("LightGBM MSE:", mean_squared_error(y, y_pred_lgb))

結果

1万行×3列のデータに対して各々のモデルで学習した結果は、自作GBDT＞sklearn＞LightGBM　の順に処理速度が速く、予測精度についてはsklearn＞自作GBDT＞LightGBM　という結果になりました。

GBDT 学習時間: 0.0023470133334437073
sklearn GBDT 学習時間: 0.0006881650001256882
LightGBM 学習時間: 0.00016488499992798703

GBDT MSE: 8363667.4633253375
sklearn MSE: 8363576.908494197
LightGBM MSE: 8364939.454835789

>>> pred
array([4918.93347788, 5004.89407361, 4966.70699032, ..., 4905.82538249,
       4905.82538249, 4905.82538249])
>>> y_pred_sklearn
array([4918.91936832, 5004.90527002, 4966.7019993 , ..., 4905.81070352,
       4905.81070352, 4905.81070352])
>>> y_pred_lgb
array([4905.08165556, 5013.22017915, 4965.05425327, ..., 4905.08165556,
       4905.08165556, 4905.08165556])
>>> y
array([3406.42980972,   74.38279199, 7711.23294004, ..., 8012.87082509,
       2804.95761046, 5334.15598591])
>>> y.mean()
4953.840769429285
>>> 

考察

これらの結果が示すことは、まず予測精度（MSE）の値から、自作GBDTのアルゴリズム構造がライブラリにかなり近い事を示しており、とくにsklearnと実質同じ計算を行っているということを数字が証明してくれています。多少の誤差が発生するのは浮動小数点の扱いによるものや、LightGBMに限っては分岐点の探索方法の違いが影響していると考えられます。
意味のない数字に対して学習を行っているため、精度としてはどのモデルも使い物にならないレベルですが、いずれのモデルも理論に忠実であると言えます。

一方、所要時間はどうしてここまでの差が出るかというと、
この差はsklearnやLightGBMに実装されているプログラミング言語や、CPU・GPUの性能などによるもので、自作のGBDTはその点作りが簡素すぎる為、データ量や計算コストに比例して時間が掛かっていると言えます。

過去にライブラリを使用せず線形回帰を実施した時に、いつまでたっても学習が終わらないという事態がありましたが、
これはどの実装を選ぶかで精度と速度が向上することを示しており、ライブラリはこれらが最適化された状態で配布されているという事ですね。
つまり、ツールの選定とはどの理論を使うかではなく、どのように内部の工夫をするのか、”理論を変えずに実務に耐えうる品質のシステムを提供する”というのが本質ではないかと考えました。

最後に、今回のモデル検証で想定していたより多くの気付きを得ることができました。

関係者の方々にお礼を申し上げます。