PySR: シンボリック回帰とは何か？

シンボリック回帰（Symbolic Regression、記号回帰とも呼ばれます）は、データを説明する数式を自動的に見つけ出す機械学習手法です。この手法では、関数の形式を事前に決めることなく、与えられたデータに最も合う数式を見つけることができます。一般的に、シンボリック回帰は遺伝的プログラミングを応用して実装されており、データ内のパターンを人間にとって理解しやすい数学的表現で捉えられるため、科学的発見やモデルの解釈性向上に役立ちます。

最近、このシンボリック回帰分野で注目を集めているのが「PySR」というオープンソースのライブラリです。PySRはPythonとJuliaを基盤にした高性能なツールで、科学的な方程式を自動的に導き出すための機能を提供します。このライブラリは、JuliaのSymbolicRegression.jlと連携して、複雑な非線形データから解釈可能な情報を抽出できる点が大きな特徴です。

今回は、シンボリック回帰の概要と、PySRの利用方法について紹介します。

シンボリック回帰の基本概念

シンボリック回帰は、データをもとに数学的な式を探索し、データのパターンを人が理解しやすい形で表現する手法です。このアプローチは、物理学や化学など、数式で自然現象を記述する必要がある分野で特に有用です。シンボリック回帰は、入力変数と出力変数の関係を表す数式を発見することで、データの背後にあるメカニズムや法則を明らかにします。

※ 通常の機械学習では、データを近似するための計算式を仮定して、その計算式のパラメータ（重み係数）を学習します。シンボリック回帰は、計算式そのものをデータから推定できます。

シンボリック回帰の仕組み

シンボリック回帰では、基本的な数学演算（加算、減算、乗算、除算、三角関数など）を組み合わせて、データに最も適合する数式を生成します。このプロセスは通常、進化的アルゴリズムや遺伝的プログラミングを用いて行われます。これらの手法を用いることで、多数の数式候補が生成され、その中から最適なものが選ばれます。シンボリック回帰の大きな利点は、得られた数式が「解釈可能」であり、関係性を単なる「ブラックボックス」モデルとしてではなく、数式として捉えられる点です。

シンボリック回帰の課題

一方で、シンボリック回帰にはいくつかの課題もあります。まず、探索空間が非常に広大であるため、計算リソースを多く消費する場合があります。また、特にノイズの多いデータでは、過剰適合（オーバーフィッティング）が生じやすく、意味のある数式を発見するのが難しくなることもあります。このため、最近の研究では、正則化手法や探索の効率を高めるアルゴリズムの工夫が進められています。

シンボリック回帰の最新技術

シンボリック回帰は近年、深層学習と組み合わせたアプローチが注目されています。例えば、ニューラルネットワークを用いて数式を生成することで、従来の手法よりも効率的かつ柔軟に数式を探索できるようになっています。また、ニューラルシンボリックモデルは、データと数式の構造を同時に学習し、より人間が理解しやすい形式で表現する試みも行われています。

シンボリック回帰の応用例

シンボリック回帰は、物理学や化学をはじめ、さまざまな分野で応用されています。たとえば、物理学では実験データから自然法則を導き出すために使われ、重力の法則や運動方程式を再発見する例もあります。化学では、化学反応の速度式を推定するために用いられ、特定の条件下での反応速度を予測することが可能です。また、経済学や生物学などでもデータ解析やモデル構築に役立てられており、特に予測モデルの理解や改善に活用されています。

PySRの紹介

PySRは、PythonとJuliaで実装されたオープンソースのシンボリック回帰ライブラリで、科学的発見を支援するために設計されています。進化的アルゴリズムを用いることで、データから解釈可能な数式モデルを効率的に発見できる高性能な分散バックエンドと柔軟な検索アルゴリズムを備えています。

このライブラリは、データの背後にある物理法則や関係を明らかにするのに役立ち、科学的データ解析やモデルの構築において広く利用されています。

PySRの基本機能

PySRは以下の基本機能を提供します：

進化的アルゴリズム: マルチ個体群進化的アルゴリズム（Multi-population Evolutionary Algorithm、MPEA）を採用し、複数の進化の経路を非同期に探索します。
ユーザー定義オペレーター: 特定の問題に合わせたオペレーターをユーザーが自由に追加できます。
分散コンピューティング: 大規模データセットへの適用も可能な分散処理機能をサポートしています。

PySRの技術的な特徴

PySRの技術的な特徴には以下が含まれます。

高性能バックエンド: JuliaのJIT（Just-In-Time）コンパイルを活用しており、計算効率が高いです。
柔軟な検索アルゴリズム: 検索アルゴリズムをカスタマイズ可能で、特定の問題に最適な数式モデルが見つけやすくなっています。

PySRのアルゴリズム

PySRでは、トーナメント選択を用いたマルチ個体群進化的アルゴリズムを採用し、突然変異や交叉によって新しい個体を生成します。進化的アルゴリズムの古典的なアプローチを改良し、より効率的な数式探索が可能です。

PySRの応用例

PySRは、さまざまな分野における科学データ解析で応用されています。

天文学: 惑星運動法則の導出などに利用されています。
材料科学: 新材料の特性予測モデルの構築に活用されています。
複雑データ解析: 他の機械学習手法と組み合わせて高度なデータ解析も可能です。

PySRの性能評価

PySRは、他のシンボリック回帰アルゴリズムや最新のディープラーニング手法と比較しても高性能であることが示されています。これは、進化的アルゴリズムを効率的に活用することで、データから解釈可能で正確な数式モデルを高速に生成できるためです。

PySRの使い方

PySRのインストールは非常に簡単で、Pythonのパッケージ管理システムであるpipを用いて行うことができます。

$ pip install pysr

では、はじめにテスト用のデータを作成します。10個の説明変数を持つデータセットをランダムに生成し、方程式 \(f(x) = x_0^2 + 2.5382 * \cos x_3 – x_7\) で目的変数を作成します。10個の変数をランダムに生成しますが、すべての説明変数を使用せずに、目的変数を説明する変数を探索できるか試したいみたいと思います。

import numpy as np
from sklearn.model_selection import train_test_split


def calc_dummy_target(x: np.ndarray) -> np.ndarray:
    return x[:, 0] ** 2 + 2.5382 * np.cos(x[:, 3]) - x[:, 7]


def load_datasets(test_size: float = 0.2) -> tuple:
    data_x = 2.0 * np.random.randn(100, 10)
    data_y = calc_dummy_target(data_x)

    return train_test_split(data_x, data_y, test_size=test_size)

次に、シンボリック回帰のインスタンスを生成する関数を定義します。PySRのシンボリック回帰用のクラスはPySRRegressorになります。多くのクラス引数があるため、ここでは説明を割愛しますが、使用ているものについてはコード中のコメントを参考にしてください。

from pysr import PySRRegressor


def build_model() -> PySRRegressor:
    return PySRRegressor(
        procs=4,  # プロセス数(=実行プロセス数): デフォルトは `cpu_count()`
        populations=15,  # 進化的プログラミングにおける集団の数: デフォルトは `15`
        population_size=33,  # 各集団における個体数: デフォルトは `33`
        ncycles_per_iteration=550,  # 各イテレーションにおける集団の10サンプルあたりに実行される突然変異の合計数:デフォルトは `550`
        niterations=5_000,  # アルゴリズムの反復回数: デフォルトは `40`
        timeout_in_seconds=60 * 1,  # 指定した秒数が経過したら探索を停止する: デフォルトは `None`
        maxsize=10,  # 方程式の最大複雑度: デフォルトは `20`
        maxdepth=10,  # 方程式の最大深度(ネスト): デフォルトは `None`

        # 利用可能な演算子の定義
        # 詳細は https://astroautomata.com/PySR/operators/ を参照のこと
        binary_operators=["*", "+", "/", "cond", "greater", "logical_or", "logical_and"],
        unary_operators=[
            "abs",
            "neg",
            "square",
            "exp",
            "cos",
            "sin",
            "tan",
            "erf",
            "inv(x) = 1/x"
        ],
        extra_sympy_mappings={"inv": lambda x: 1 / (x + 1.0E-8)},

        complexity_of_constants=0.5,  # 定数に対する複雑度を指定: デフォルトは `1`
        select_k_features=4,  # 使用する特徴量の数を指定: デフォルトは `None`
        progress=True,  # 進行状況バーを使用するかどうか: デフォルトは `True`
        weight_randomize=0.00023,  # 突然変異による方程式破棄で方程式が新たにランダムに生成される相対的な可能性: デフォルトは `0.00023`
        cluster_manager=None,  # 分散コンピューティングの設定 ("slurm", "pbs", "lsf", "sge", "qrsh", "scyld", or "htc"): デフォルトは `None`
        precision=32,  # 計算精度 (16 or 32 or 64): デフォルトは `32`
        warm_start=False,  # `fit`が繰り返し呼ばれたときに前回の結果から再開するかどうか: デフォルトは `False`
        turbo=True,  # (Experimental) LoopVectorization.jl を使用して検索評価を高速化するかどうか: デフォルトは `False`
    )

必要な定義が準備できたので、上記の関数を利用して実際にシンボリック回帰を実行します。最後に精度を確認していますが、指標はR²スコアです。結果を見ると、訓練データ、テストデータともに非常に高い精度になっていることが確認できます。

# データセットの読み込み
train_x, test_x, train_y, test_y = load_datasets()

# モデルを作成
model = build_model()

# 学習
model.fit(train_x, train_y)

# 精度
print("Train score: ", model.score(train_x, train_y))
# => Train score:  0.9999999994808751
print("Test score: ", model.score(test_x, test_y))
# => Test score:  0.9999999993663083

次に、どのような方程式を使ったのか確認します。

print(model)

モデルの詳細を表示すると、次のようになりました。

PySRRegressor.equations_ = [
           pick      score                                         equation          loss  complexity
        0         0.000000                                               x3  6.398231e+01           1
        1         2.571436                                       square(x0)  4.889893e+00           2
        2         0.015624                     abs(square(x0 + 0.09493776))  4.739453e+00           4
        3         0.435242                             square(x0) + neg(x7)  3.066937e+00           5
        4         0.077010               (neg(x7) * 0.7502878) + square(x0)  2.839618e+00           6
        5         0.034463   square(x0) + ((-0.41587457 + x7) / -1.3405459)  2.743423e+00           7
        6         0.890212                 (square(x0) + cos(x3)) + neg(x7)  1.126364e+00           8
        7         0.532187            (square(x0) + tan(cos(x3))) + neg(x7)  6.615353e-01           9
        8  >>>>  17.116185  (square(x0) + (cos(x3) / 0.39401513)) + neg(x7)  2.438307e-08          10
]

LaTeX形式でも表示して、想定している方程式を探索できているか見ていきます。すると、探索した方程式は、\(x_{0}^{2} – x_{7} + \frac{\cos{\left(x_{3} \right)}}{0.394}\)となっていることが分かります。\(\frac{1}{0.394} \approx 2.538\) であるため、ランダムに作成したデータセットの方程式と一致しています。

print(model.latex())
# => x_{0}^{2} - x_{7} + \frac{\cos{\left(x_{3} \right)}}{0.394}

まとめ

今回は、シンボリック回帰とPySRを紹介し、データから解釈可能な数式を導く方法を解説しました。シンボリック回帰は、データサイエンスにおいて変数間の関係を理解するための強力なツールです。

またシンボリック回帰のライブラリであるPySRは、シンボリック回帰の仕組みを自動化し、解釈可能な予測モデルを導き出すのを支援します。今後、PySRとディープラーニングの統合が期待され、さらに多様なデータセットへの適用が見込まれています。PySRの詳細はGitHubリポジトリやAPIリファレンスを参照してください。

More Infomation:
arXiv:2305.01582, Miles Cranmer (Princeton University and Flatiron Institute), 「Interpretable Machine Learning for Science with PySR and SymbolicRegression.jl」, https://arxiv.org/abs/2305.01582

Blog